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folecular  dynamics  are  computed  for  a  model  S*2  reaction 
CfTj+  CrijCl  -*  Cioij  +  CfHjin  water  and  are  found  to  be  strongly  dependent  on  the 
instantaneous  local  configuration  of  the  solvent  at  the  transition  state  barrier.  There  are 
significant  deviations  from  the  simple  picture  of  passage  over  a  free  energy  barrier  in 
the  reaction  coordinate,  and  thus,  a  marked  departure  from  Transition  State  Theory 
occurs  in  the  form  of  barrier  recrossings.  Factors  controlling  the  dynamics  are  dis¬ 
cussed.  and,  in  particular,  die  rate  of  change  of  atomic  charge  distribution  along  the 
reaction  coordinate  is  found  to  have  a  major  effect  on  the  dynamics.  A  simple  frozen 
solvent  theory  involving  nonadiabatic  solvation  is  presented  which  can  predict  the  out¬ 
come  of  a  particular  reaction  trajectory  by  considering  only  the  interaction  with  the 
solvent  of  the  reaction  system  at  the  gas-phase  transition  barrier.  The  frozen  solvent 
theory  predicts  the  adjustment  (the  transmission  coefficient  k)  needed  to  make  the 
transition  state  theory  rate  agree  with  the  outcome  of  the  molecular  dynamics  trajec- 
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I.  Introduction 

The  bimolecular  nucleophilic  substitution  SN2  reaction,  X~  +  RY  ->  RX  +  Y~,  is  one  of  the  most 
fundamental  of  organic  reactions.  It  has  been  extensively  studied  in  solution,  and  these  studies  are  a 
source  of  many  basic  ideas  in  physical-organic  chemistry;  examples  include  nucleophilicity,  leaving- 
group  ability  and  solvent  effects.1'3  This  reaction  is  well  known  to  be  experimentally  sensitive  to  sol¬ 
vent  polarity,1'4  and,  beginning  with  the  famous  work  of  Hughes  and  Ingold,3'7  solvent  effects  have 
been  interpreted  in  terms  of  the  equilibrium  solvation  thermodynamic  free  energies  of  activation  of  the 
transition  state  complex.  This  theme  has  been  pursued  in  modem  statistical  mechanics  in  a  Monte  Carlo 
equilibrium  simulation  of  Cl~  +  CH3C1  by  Jorgensen  and  coworkers8-9  and  in  integral  equation  work  by 
Chiles  and  Rossky.10  In  this  paper  we  use  the  Molecular  Dynamics  simulation  technique  to  explore  a 
quite  different  and  novel  aspect  of  the  reaction  pathway:  The  role  of  polar-solvent  dynamics  and 
configurations  in  affecting,  at  the  molecular  level,  the  reaction  trajectories  and  the  value  of  the  rate  con¬ 
stant. 

These  last  ingredients  have  been  explored  for  general  charge  transfer  reactions  in  analytic  model 
studies  by  van  der  Zwan  and  Hynes.11'13  These  studies  point  strongly  to  the  importance  of  nonequili¬ 
brium  solvation  effects  on  charge  transfer  rates  in  polar  solvents.  The  present  work  examines  these 
questions  for  a  more  realistic  reaction  in  a  realistic  solvent. 

Our  method  is  simple  and  direct  and  follows  along  the  general  lines  of  our  previous  investiga¬ 
tions14- 13  of  the  A  +  BC  reaction  in  rare  gas  solvent.  We  use  molecular  dynamics  computational  simu¬ 
lation  to  examine  model  SN2  reactions  roughly  patterned  after  the  Cl~  +  CH3C1  reaction  in  water.  The 
intrinsic  SN2  potential  is  taken  to  be  the  form  of  a  double  well,  established  in  gas  phase  experi¬ 
ments2- 3-  16  and  quantum  mechanical  calculations,8-17'20  but  we  take  the  methyl  group  as  simplified  to  a 
single  united  atom.  A  fairly  realistic  model  for  flexible  water  is  adopted.21 

We  find  that,  even  for  high,  sharp  central  barriers,  the  ability  of  the  solvent  to  influence  the  tra¬ 
jectories  of  the  reacting  atoms,  and  thus  to  lower  the  rate  constant  from  the  simple  Transition  State 
Theory  prediction,  is  significant.  The  molecular  level  picture  of  the  reaction  dynamics  differs  consider¬ 
ably  from  that  painted  by  an  equilibrium  free  energy  picture  in  which  equilibrium  solvation  is  assumed. 
We  find  that  the  transition  state  is  characterized  by  nonequilibrium,  nonadiabatic  solvation;11'13  the  sol¬ 
vent  is  effectively  "frozen"  and  cannot  adjust  rapidly  enough  on  the  relevant  time  scales  to  provide 
equilibrium  solvation.  We  explore  this  reaction  as  a  function  of  reaction  barrier  height  and  curvature, 
characteristics  of  charge  redistribution  along  the  reaction  coordinate,  and  factors  controlling  solvent 
molecule  reorientation  time. 
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The  outline  of  this  paper  is  as  follows.  In  Sec.  II  we  describe  the  methodology,  selection  of  die 
potential  energy  surfaces  and  the  prescription  for  the  charge  redistribution  rate  along  the  reaction  coor¬ 
dinate.  In  Sec.  m  we  analyze  the  reaction  results  for  a  barrier  height  similar  to  that  of  the  actual 
CT  +  CHjCl  reaction,  while  Sec.  IV  explores  the  consequences  of  different  barrier  heights  and  curva¬ 
tures,  variation  of  charge  switching  rate  and  modification  of  the  water  hydrogen  mass.  Section  V  is  the 
conclusion  and  summary.  Certain  details  of  the  calculation  are  relegated  to  the  Appendices.  A  com¬ 
panion  paper22  describes  in  more  detail  an  analytic  model  to  describe  our  results. 

n.  Methodology  and  Potentials 


A.  Molecular  Dynamics 


Molecular  dynamics  is  used  to  compute  time  histories  of  the  position  and  momentum  coordinates 
for  all  of  the  atoms  in  a  system  subject  to  truncated  octahedral23  periodic  boundary  conditions.  This  is 
accomplished  by  numerical  integration  of  Hamilton’s  equations  of  motion  for  classical  particles 


(2.1) 


and 


P* 


a h 

'  dr*  ’ 


(2.2) 


in  which  pN  and  r"  are  the  conjugate  momentum  and  position  variables  for  a  system  containing  N 
unique  particles,  and  H  is  the  full  Hamiltonian  for  the  system. 

A  modified  Verlet  algorithm,24"26  incorporating  a  time  step  of  0.5  femtoseconds,  is  used  to 
integrate  these  equations  of  motion.  There  are  64  water  molecules  used  in  most  of  the  liquid  simula¬ 
tions,  and  263  water  molecules  in  a  limited  set  of  test  calculations. 


These  calculations  are  carried  out  on  a  Floating  Point  Systems  AP120B  array  processor  attached 
to  a  VAX  11/780  host  computer.  For  64  water  solvent  molecules  and  the  reacting  solute,  1000  time 
steps  require  300  seconds  of  real  time  on  the  array  processor.  The  calculations  presented  in  this  paper 
required  80  days  of  array  processor  time,  and  would  have  required  7.S  years  of  dedicated  VAX  11/780 
time  without  the  array  processor. 


B.  Potentials 

The  potential  energy  surface  used  in  the  Hamiltonian  in  Eqs.  (2.1)  and  (2.2)  is  intended  to  be  a 
simplified  representation  of  the  model  Sn2  reaction,  Cl"  +  CH3C1  ->  CICH3  +  Cl",  occurring  in  liquid 
water  solvent  The  total  potential  is  composed  of  the  potentials  representing  smaller  parts  of  the  prob¬ 
lem  as  described  below. 


1.  S*2  Gas  Phase  Potential 

Our  Sn2  gas  phase  potential  is  very  roughly  modeled  after  the  one-dimensional 
G~  +  CH3CI  — »  CHjCl  +  Cl“  ab  initio  gas  phase  potential  of  Jorgensen  and  co-workers.8*9  The  Jorgen¬ 
sen  potential  was  only  given  along  a  reaction  coordinate  defined  as  the  gas  phase  minimum  potential 
energy  path  between  reactants  and  products.  To  describe  a  full  three-dimensional  potential  energy  sur¬ 
face  for  our  model  SN2  reaction  we  use  a  London-Eyring-Polanyi-Sato  (LEPS)27  potential  energy  func¬ 
tion.  This  incorporates  the  main  features  of  the  one-dimensional,  double-well  potential  determined  by 
Jorgensen,  including  the  calculated  barrier  height,  but  uses  a  united  atom  representation  for  the  central 
CHj  group.  The  form  of  the  LEPS  potential  is 

*“"0Ws)  =  fit  +  Qi  +  Cj  -  (A  +  A  +  A  -  Vz  -  V i  -  Vi)*  .  (2.3) 

in  which  ri  and  r2  are  the  variable  bond  lengths  between  the  CHj  united  atom  and  the  two  G  atoms,  r3 
is  the  bond  length  between  the  two  G  atoms,  and  Q,  and  J,  are  linear  combinations  of  the  "singlet"  and 
"triplet"  diatomic  energies: 

Q,<r,) = ('£,  +  %n , 
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J,<r.)  =  (%  ~  %)f2  •  (2.5) 

Hoe,  the  'singlet*  ground  state  energy  is  represented  as  a  Morse  potential, 

'EM  =  >a|i  -  exp  pft(r4  -  -  lD, ,  (2.6) 

while  the  "triplet"  diatomic  energy  is  an  "anti-Morse"  potential, 

'EM  =  *a|i  +  exp[-sp((r,  -  Jrf)]|  -  JD, .  (2.7) 

The  LEPS  potential  reduces  to  the  diatomic  bound  Morse  potential  for  an  infinitely  separated 
atom  and  diatomic  ( i.e .  CHj— Cl)  and  is  highly  adaptable,  containing  18  non-symmetry  adapted  param¬ 
eters  (‘Dj,  J£>j,  'p,,  sp„  ‘if,  Jif,  for  «  «  1,2,3)  which  allow  the  simultaneous  adjustment  of  barrier  height 
and  barrier  frequency,  to*  defined  as  the  absolute  value  of  the  imaginary  frequency  of  an  inverted  para¬ 
bolic  approximation  to  the  potential  at  the  saddle  point  along  the  minimum  energy  pathway  from  reac¬ 
tants  to  products,  i.e.  die  reaction  coordinate.  We  report  all  frequencies  in  spectroscopic  wavenumber 
units,  i.e.  (cy2rtc)  cm'1,  to  facilitate  comparison  with  spectroscopic  vibrational  properties.  Table  1 
shows  the  values  of  the  18  LEPS  parameters  used  in  the  simulation  of  our  model  SN2  reaction  and  vari¬ 
ations  which  incorporate  the  basic  features  of  the  double-well  potential  employed  by  Jorgensen8  and  the 
isolated  Morse  potential  fitted  to  die  equilibrium  value28  of  the  normal  mode  asymmetric  stretch  vibra¬ 
tional  frequency,  732.8  cm-1,  of  CHj — Cl.  (Note  that  the  LEPS  potential  per  se  allows  bonding 
between  the  two  chlorine  atoms.  We  choose  the  Cl-Cl  parameters  to  be  free  variables  for  fitting  our 
potential  to  Joigensen’s  potential,  so  they  do  not  represent  physical  parameters.)  Figure  1  illustrates 
potential  energy  contours  for  the  13.9  kcal/mol  barrier  in  which  die  axes  are  skewed  44.4°  from  perpen¬ 
dicular  to  define  a  coordinate  system  which  diagonalizes  the  kinetic  energy  matrix  for  this  system.27-29 
(This  allows  classical  dynamics  to  be  represented  by  a  single  point  mass  sliding  on  the  potential  energy 
surface.)  The  key  features  of  this  potential  are  two  wells  characterizing  ion-molecule  complexes 
separated  by  a  high  transition  barrier.  Most  of  our  calculations  are  performed  on  this  three  dimensional 
LEPS  potential,  which  has  a  barrier  of  13.9  kcal/mol  with  respect  to  the  bottom  of  the  ion-dipole  com¬ 
plex  wells  and  3.6  kcal/mol  with  respect  to  separated  CICHj  +  Cl~. 

2.  Solvent  Potential 

We  adopt  the  semiempirical  flexible  water  intermolecular  potential  of  Watts21  combined  with  the 
spectroscopic  intramolecular  water  potential  of  Kuchitsu  and  Morino.21-30  These  water  potentials  have 
been  used  in  molecular  dynamics  simulations  to  calculate  thermodynamic  quantities,  energy,  free  energy 
and  constant  volume  heat  capacity,  and  vibrational  normal  mode  frequencies  31-32  in  good  agreement 
with  experiment  after  a  quantum  correction  is  added.  In  addition,  the  rotational  time  correlation  func¬ 
tions  and  radial  distribution  functions  of  these  waters  compare  favorably  with  other  water  models  33 

3.  Solvent-Solute  Potential 

Charges  on  the  reaction  system  interact  with  charges  on  water  molecules  through  coulomb  poten¬ 
tials.  The  dipole  moment  on  water  (-1.86  D)  is  set  by  assigning  partial  charges  of  -0.66  electrons  to 
oxygen  atoms  and  +0.33  electrons  to  hydrogen  atoms.21  Because  the  waters  are  internally  flexible,  the 
dipole  moment  can  fluctuate  with  changing  internal  nuclear  coordinates.  The  partial  charges  on  the 
dune  atoms  in  our  model  SN2  reaction  system  sum  to  give  the  overall  negative  charge  of  one  electron. 
The  distribution  of  this  charge  along  the  reaction  coordinate  varies  smoothly  from  reactants  where  the 
net  charge  is  on  the  attacking  chloride  ion,  to  the  transition  barrier  where  it  is  shared  equally  between 
the  two  chlorine  atoms,  to  the  products  where  it  is  localized  on  the  leaving  chloride.  Figure  2  illustrates 
the  distribution  of  the  partial  charges  on  the  three  atoms  as  a  function  of  a  conveniently  defined  (for 
computational  purposes)  reaction  measure,  rjj  -  rjc,  in  which  is  the  distance  between  one  G — CH3 
atom  pair  and  rtc  is  the  distance  between  the  other  CH3 — Cl  atom  pair.  This  variation  is  largely 
modeled  after  the  work  of  Jorgensen,8  but  with  some  modifications  of  the  charge  switching  parameters 
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Tabic  1.  [Cl — CHj — Cl]-  LEPS  Potential  Parameters 


Barrier  Height 
(with  respect  to  the 
bottom  of  the  well, 
kcal/mol) 

31.9 

13.9 

4.9 

Barrier 
Frequency 
(magnitude,  cm-1) 

894 

461 

232 

'r?  (A)  i  -  U 

1.776382 

1.776382 

1.776382 

M  <*> 

2.094857 

2.094857 

2.094857 

M  (k)  i  -  u 

1.789255 

1.785014 

1.795583 

>rS  <*> 

2.135243 

2.186060 

2.175774 

%  (A"1)  i  -  1,2 

0.934381 

0.929968 

0.936038 

%  (A-1) 

0.639340 

0.432955 

0.705479 

3P,  (A-1) «  -  1,2 

4.512772 

4.822681 

4.853273 

’fh  (*-■> 

1.048338 

1.016811 

1.036805 

lDt  (kcal/mol)  t  -  1,2 

233.134491 

234.524674 

232.308932 

iDi  (kcal/mol) 

308.573800 

64.925971 

-88.274115 

SD,  (kcal/mol)  <  -  1,2 

261.994490 

220.244820 

261.021681 

3D3  (kcal/mol) 

274.953525 

284.999867 

274.257166 

for  the  central  united  atom  (cf  Appendix  A).  This  shifting  charge  on  the  solute  interacts  with  the  sur¬ 
rounding  solvent  molecules  through  Coulombic  potentials,  with  the  partial  charge  on  each  solvent  atom 
fixed. 

A  moving  charge  in  the  reacting  system  has  an  interesting  effect  on  the  forces  between  the  solute 
and  solvent  atoms.  The  three  atoms  in  the  reacting  system  will  "feel"  a  force  from  the  solvent  related  to 
the  switching  charge  which  varies  along  the  reaction  coordinate  in  addition  to  the  normal  Coulombic 
force  for  fixed  charges.  Near  the  transition  barrier  the  charge  is  shifting  rapidly  between  the  attacking 
and  leaving  chloride  ions,  and  the  extra  force  which  arises  is  proportional  to  the  rate  of  change  of  this 
charge  distribution  with  respect  to  the  reaction  coordinate.  At  any  instant  only  the  reaction  system  feels 
this  force,  while  all  solvent  molecules  are  unaffected  by  this  extra  charge  switching  force;  the  total 
force  between  the  solute  and  solvent  molecules  is  nevertheless  conservative.  This  force  is  a  "heavy  par¬ 
ticle  polarization  force",  since  displacement  in  the  reaction  coordinate  shifts  the  electronic  distribution 
instantaneously  (in  die  electronic  adiabatic  limit).  We  will  find  that  this  shifting  charge  force  has 
important  consequences  in  the  dynamics  of  the  SN2  reaction. 

The  van  der  Waals  forces  between  solute  and  solvent  atoms  are  described  using  simple  Lennard- 
Jones  6-12  potentials.  If  required,  the  Lennard-Jones  parameters  for  the  solute  could  be  made  to  vary 
with  reaction  coordinate*  in  a  similar  manner  to  the  charge  variation,  but  we  choose  to  keep  them  con¬ 
stant  Mathematical  details  of  die  solvent-solute  potential  and  forces  are  given  in  Appendix  A. 

The  computed  G-H  and  G-O  radial  distribution  functions  (rdf  s)  of  the  reactants  at  the  transition 
barrier  with  the  water  solvent  compare  well  to  the  rdf  s  of  Jorgensen.9  In  addition,  the  pair  interaction 
distribution  functions  for  a  single  chloride  ion  in  water,  computed  with  the  pure  water  potential  (see 
above)  and  with  the  Lennard-Jones  and  Coulombic  parameters  for  chlorine-water  interactions,  compare 
well  with  die  pair-interaction  energies  determined  by  Vogel  and  Heinzinger  for  lithium  chloride  and 


cesium  chloride  in  water.34’35 


4.  Switching  Functions 

In  simulating  this  Sn2  reaction,  CT  +  CHjCl  -*  C1CH3  +  CP  in  a  polar  solvent,  special  care  must 
be  taken  in  ’feathering”,  i.e.  smoothly  ramping,  the  potentials  properly  to  zero  by  the  edge  of  the  boun¬ 
dary  conditions  to  avoid  force  discontinuities  and  nonphysical  buildup  of  net  charges  on  molecules.  To 
accomplish  this  we  define  certain  switching  regions  wherein  the  forces  are  damped  near  the  edge  of  our 
periodic  boundary  conditions  on  a  molecule-by-molecule  basis.  The  interactions  among  solvent 
molecules  and  of  the  solvent  molecules  with  the  charged  reactant  species  are  computed  for  each 
molecule-molecule  pair.  The  distance  between  a  particular  molecule-molecule  pair  is  defined  as  die  dis¬ 
tance  between  one  molecule  and  the  nearest  periodic  image  of  the  other  molecule23  (cf  Appendix  B). 
Forces  calculated  between  these  two  imaged  molecules  are  damped  when  this  distance  becomes  close  to 
the  periodic  boundary  length,  and  are  zero  when  this  distance  is  equal  to  the  periodic  boundary  length. 

A  difficulty  in  simulating  the  motion  of  the  solute  atoms  along  the  reaction  coordinate  lies  in 
maintaining  the  overall  neutrality  of  the  methyl  chloride  species  away  from  the  transition  barrier.  On 
the  reactant  side,  (O'  +  CH3C1),  die  CH3Q  is  imaged  with  all  other  molecules  as  a  single  molecule 
because  it  has  a  non-zero  dipole  moment.  Similarly,  on  the  product  side  (OCH3  +  O-),  the  ClCHj  is 
imaged  as  one  molecule.  Although  imaging  the  entire  [O— CH3 — Cl]“  complex  is  possible,  the  short 
range  interactions  with  either  the  ion  or  the  molecule  would  be  underestimated  due  to  the  damping  by 
the  switching  function  at  large  distances.  To  solve  this  problem  the  ion  and  the  molecule  are  imaged 
separately  by  switching  the  definition  of  the  grouping  from  reactant  to  products  on  either  side  of  the 
transition  barrier  in  a  continuous  manner.  This  is  explained  in  detail  in  Appendix  B. 


C.  Initial  Conditions 


The  molecular  dynamics  of  activated  barrier  crossing  are  computed  using  a  technique  due  to 
Keck36  and  Anderson,37  (see  also  Bennett3*  )  in  which  trajectories  are  initialized  from  an  equilibrium 
distribution  at  the  top  of  the  energy  barrier,  and  dynamics  are  computed  both  forward  and  backward  in 
time  from  the  initial  conditions.  This  general  method  has  been  applied  to  the  study  of  defect  motion,39 
isomerization  dynamics,40  protein  dynamics,4M3  and  surface  desorption,44  as  well  as  to  the  dynamics 
of  the  A  +  BC  reaction  in  rare  gas  solvents.14’ 15 


We  use  molecular  dynamics  to  compute  an  ensemble  of  trajectories  initially  selected  at  the  transi¬ 
tion  barrier  but,  as  described  below,  originating  and  ending  with  reactants  or  products  and  analyze  the 
results  by  examining  statistical  and  average  properties  of  the  ensemble.  Examples  of  these  are  the 
number  of  barrier  recrossings  and  the  details  of  the  energy  exchange  with  the  solvent 

The  average  value  of  some  quantity,  (A),  in  a  classical  ensemble  is 


M)  = 


JJ 


(2.8) 


where  r*  is  a  vector  of  the  N  Cartesian  coordinates  of  the  system,  pN  is  a  vector  of  the  momentum 
coordinates  conjugate  to  r*.  H{ r*,p*)  is  tire  full  Hamiltonian  of  the  system,  kB  is  Boltzmann’s  constant 
and  T  is  the  temperature.  The  initial  conditions  r*  and  p*  are  selected  from  a  Boltzmann  distribution 
on  the  transition  barrier  which,  in  the  gas  phase  (and  for  very  weak  solvent-solute  interactions),  can  be 
conveniently  defined  as  lying  on  the  symmetric  stretch  and  bend  coordinates  of  [Cl — CH} — Cl]' 
through  the  lowest  lying  reddle  point  on  the  gas  phase  LEPS  potential  surface.  The  reaction  coordinate 
is  defined  to  be  the  asymmetric  stretch,  and  the  absolute  value  of  the  normal  mode  imaginary  frequency 
associated  with  the  asymmetric  stretch  coordinate  on  the  gas  phase  potential  surface  is  defined  as  the 
barrier  frequency,  to*.  Although  we  describe  in  a  previous  paper15  a  method  for  evaluating  Eq.  (2.8) 
under  conditions  of  weak  solvent-solute  interactions,  that  method  is  not  applicable  here  since  there  are 
very  strong  charge-charge  interactions  between  the  solute  and  water  solvent 


To  describe  our  current  procedure,  we  first  note  that  the  reaction  is  symmetric  with  respect  to 
reactants  and  products.  This  implies  that  the  average  value  of  the  asymmetric  stretch  coordinate  is  zero. 
Thus,  if  there  is  a  single  saddle  point  the  average  location  of  the  transition  barrier  must  correspond  to 
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zero  asymmetric  stretch  and  lie  somewhere  along  the  symmetric  stretch  and  bend  coordinates  of 
[Cl — CHj — Cl]~.  The  potential  energy  surface  on  which  we  choose  initial  conditions  for  the  solute 
must  include,  in  addition  to  the  gas  phase  potential,  the  effect  of  the  strong  intermolecular  potential 
energy  of  interaction  with  die  solvent  Therefore,  our  initial  conditions  are  chosen  by  first  constraining 
the  asymmetric  stretch  position  and  momentum  of  [Cl — CH3 — Cl]~  to  be  zero  while  molecular  dynam¬ 
ics  selects  values  of  all  the  other  degrees  of  freedom  in  the  system.  This  constraint  is  performed  by 
forcing,  at  each  time  step,  the  lengths  of  the  tvo  CH3 — Cl  bonds  to  be  equal,  and  removing  all 
momenta  and  force  along  the  asymmetric  stretch.  After  a  2.0  ps  equilibration  period  to  allow  the  other 
degrees  of  freedom  to  randomize,  the  asymmetric  stretch  momentum  is  initialized  from  a  Boltzmann 
distribution  at  temperature  298  K,  the  constraint  on  the  reaction  coordinate  is  released,  and  dynamics 
are  computed  both  forward  and  backward  in  time  long  enough  to  observe  on  both  ends  of  the  trajectory 
whether  a  stable  product  or  reactant  is  formed.  It  is  important  to  note  that,  in  contrast  to  our  earlier  A 
+  BC  in  rare  gas  solvent  studies14-15  where  we  used  quadrature  techniques  which  gave  each  trajectory  a 
different  weighting  in  the  ensemble,  in  this  paper  we  let  the  dynamics  randomly  choose  the  initial  con¬ 
ditions  so  that  all  of  the  trajectories  have  equal  weighting  in  the  ensemble  average. 

To  investigate  the  possibility  that  die  solvent  might  be  "trapped"  in  a  local  solvent  potential 
minimum  (and  thus  that  a  similar  solvent  configuration  might  persist  between  trajectory  runs),  a  series 
of  tests  was  made  in  which  the  entire  [Cl — CH3 — Cl]*  molecule  is  held  fixed  while  the  solvent  is 
heated  to  a  temperature  of  5000  K,  the  molecular  dynamics  run  for  1.0  ps  of  simulation  time  in  order  to 
expel  the  solvent  out  of  any  local  minima,  and  the  solvent  then  cooled  back  down  to  298  K.  The  high 
temperature  equilibration  is  found  to  give  results  similar  to  those  obtained  with  298  K  equilibration,  and 
the  majority  of  our  runs  are  thus  made  with  equilibration  at  298  K. 

m.  Results  for  Model  Cl~  +  CH3C1  Reaction 

In  this  section  we  present  the  results  for  the  model  SN2  reaction  of  Cl*  +  CH3C1  in  both  gas  phase 
and  in  liquid  water.  There  are  64  water  molecules  in  the  simulation  of  the  liquid  at  a  density  of  1.0 
g  cm*3.  Using  a  larger  system  is  desirable,  but  this  would  make  the  calculations  unpractically  time- 
consuming  since  the  amount  of  computer  time  required  to  simulate  the  dynamics  is  considerable,  even 
for  the  small  solvent  system  of  64  waters.  Thus,  a  small  subset  of  the  dynamics  is  computed  for  a 
larger  system  containing  263  solvent  water  molecules  to  check  that  there  are  no  significant  differences 
between  the  dynamics  obtained  for  the  smaller  and  larger  systems.  Molecular  dynamics  calculations 
with  the  larger  263  water  system  take  approximately  16  times  longer  than  with  the  smaller  system,  so  it 
is  possible  to  generate  only  6  trajectories  (with  a  total  simulated  time  of  15  ps)  per  day,  even  with  the 
use  of  the  array  processor. 

A.  Gas  Phase  Dynamics 

Although  our  focus  is  on  reaction  dynamics  in  solution,  the  dynamics  in  the  gas-phase2-45  is  of 
sufficient  interest  to  warrant  a  brief  exposition.  Figure  1  shows  an  energy  contour  plot  of  the  gas-phase 
potential  energy  surface  in  which  the  double  well  character  centered  about  the  asymmetric  stretch  reac¬ 
tion  coordinate  in  the  transition  barrier  region  can  be  clearly  seen.  The  central  transition  energy  barrier 
is  higher  in  energy  than  the  asymptotic  potential  energy  of  the  separated  Cl*  ion  and  CH3C1  by  3.6 
kcal/mol,  so  there  is  sufficient  energy  for  a  trajectory  starting  at  the  transition  barrier  to  escape  from  the 
low  energy  product  well  past  the  transition  barrier  and  form  separated  products.  However,  this  is  not 
what  is  observed;  instead  the  trajectories  are  trapped  in  the  ion-dipole  complex  energy  well  for  a  time 
measuring  much  longer  than  our  standard  simulation  time  of  ±  0.5  ps.  A  vocationally  "hot"  complex  is 
formed  which  oscillates  inside  the  well  for  some  time  before  breaking  up.  In  addition,  there  are  some 
recrossings  of  the  transition  barrier  as  the  complex  passes  from  the  product  well  to  the  reactant  well. 
Table  2  shows  crossing  patterns  for  128  trajectories,  each  of  which  is  integrated  out  to  time  0.5  ps 
before  and  after  the  transition  barrier.  Here,  "products”  are  defined  as  the  C1CH3 — Cl*  ion-molecule 
complex.  Trajectories  always  begin  with  initially  forward  momentum  toward  products  along  the  asym¬ 
metric  stretch  reaction  coordinate,  and  Table  2  indicates  that  107  out  of  128  trajectories  became  "pro¬ 
ducts"  on  die  ±  0.5  ps  time  scale  on  which  we  followed  the  dynamics.  However,  very  few  of  these  tra¬ 
jectories  actually  escaped  the  wells,  but  instead  formed  a  metastable  complex.  The  other  21  trajectories 
actually  recossed  the  transition  barrier  at  least  once  and  then  remained  localized  in  the  appropriate 


Number  of  Crossings 

Initial  and  Final  States  of  Trajectories 

Reactant  -  Reactant 

Reactant  -  Product 

Product  -  Reactant 

Product  -  Product 

1 

0 

107 

0 

0 

2 

10 

0 

0 

9 

3 

0 

0 

2 

0 

complex  well. 

One  interesting  way  to  examine  the  features  of  isolated-complex  dynamics  is  to  monitor  the  redis¬ 
tribution  of  internal  energy  in  the  reaction  system  as  a  function  of  time.  One  way  to  partition  the  total 
energy  is  into  [Cl — CH3— -Cl]“  potential  energy,  [Cl — CH3 — Cl]"  potential  energy  plus  CH3 — Cl  vibra¬ 
tional  lrinetic  energy,  CH3 — G  rotational  and  translational  kinetic  energy,  and  the  translational  kinetic 
energy  of  Cl".  Figure  3  shows  average  energy  partitioning  patterns  as  a  function  of  time  for  128  trajec¬ 
tories  in  the  gas  phase.  The  vibrational  kinetic  energy  plus  total  [Cl — CH} — Cl]"  potential  energy  is 
seen  to  increase  to  a  maximum  in  about  0.2S  ps  before  and  after  the  transition  barrier  but  decreases  at 
longer  times.  This  build  up  in  energy  and  subsequent  decrease  is  observed  on  a  shorter  time  scale  for 
the  translational  and  rotational  degrees  of  freedom.  The  kinetic  energy  at  longer  times  is  seen  to 
transform  into  potential  energy,  and  a  metastable  complex  is  formed.  Complex  formation  is  a  result  of 
the  strong  charge-dipole  interaction  close  to  the  transition  barrier  and  is  characterized  by  die  double¬ 
well  feature  of  the  potential  surface.  Even  though  there  is  total  energy  in  excess  of  that  required  to 
escape  to  separated  products  from  the  bottom  of  the  ion-dipole  well,  there  clearly  are  coupled  non¬ 
reactive  degrees  of  freedom  in  the  system  into  which  the  energy  is  redistributed  on  the  time  scale  of  our 
calculation.  (Note  that  the  methyl  umbrella  vibration  is  not  a  dynamical  variable  here.) 

While  the  isolated-complex  dynamics  is  an  interesting  subject  for  future  research  with  implica¬ 
tions  for  RRKM  theory,  we  present  it  here  solely  to  contrast  with  the  reaction  in  solution.  The  follow¬ 
ing  section  examines  the  dynamics  of  this  reaction  system  in  water  solvent. 

B.  Dynamics  in  Aqueous  Solution 

In  the  standard  view5-6  of  solvent  effects  on  the  rates  of  SN2  reactions,  the  solvent  is  assumed11' 

13  to  be  in  equilibrium  at  each  stage  of  the  "intrinsic”  X"  +  RY  reaction  coordinate.  In  fact,  it  is 
assumed  that  the  reaction  coordinate  depends  only  on  the  X"  +  RY  system,  and  the  polar  solvent  only 
modifies  the  gas-phase  free  energy  profile  along  the  reaction  coordinate  by  lowering  the  free  energy  of 
the  higher-charge-density  reactants  X”  +  RY  and  products  RX  +  Y"  more  than  the  free  energy  of  the 
dispersed-charge-density  transition  state  X  *" — R  ** — Y  Thus,  the  free  energy  of  activation  increases 
in  solution  relative  to  the  gas  phase, and  die  deep  wells  on  either  side  of  the  transition  barrier  diminish 
considerably.2* 3,8,9  This  is  illustrated  in  Fig.  4,  which  shows  a  schematic  representation  of  the  free 
energy  profile  along  the  reaction  coordinate  for  a  symmetric  SN2  reaction  of  the  form  X"  +  RY  in  both 
the  gas  phase  (solid  line)  and  in  a  polar  solvent  (dashed  line).2,3,5,6  These  features  are  borne  out  in 
explicit  statistical  mechanical  calculations  of  reaction  free  energies.*'10 

The  free  energy  picture  is  a  restricted  view  which  makes  no  allowances  for  possible  dynamical  or 
non-equilibrium  effects  of  solvent  molecules  on  the  details,  or  even  the  fate  of  individual  reaction  tra¬ 
jectories.  Does  die  equilibrium  picture  for  the  solvent  hold  at  all  points  along  the  reaction  coordinate, 
and  does  the  reaction  always  proceed  to  formation  of  products  upon  reaching  the  transition  barrier  with 
positive  flux?  Equivalently,  are  there  nonequilibrium  solvation  effects11'13  which  lead  to  recrossing  of 
the  TST  surface?  The  answers  to  these  questions  are  critical  in  understanding  SN2  reaction  dynamics  in 
solution. 
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1.  Barrier  Recrossings 

Each  trajectory  has  initially  positive  momentum  in  the  asymmetric  stretch  reaction  coordinate 
which  sends  it  toward  products  GCHj  +  Cl-.  We  have  examined  400  solution  phase  trajectories  with 
with  a  solvent  composed  of  64  water  molecules.  We  find  that  the  "no  recrossing"  assumption  is  often 
violated;  further,  the  reaction  fate  is  typically  decided  within  the  very  short  time  period  of  0.02  ps.  The 
13.9  ltcal/mol  energy  barrier  crossing  patterns  are  presented  in  Table  3  and  show  that  trajectories  in 
solution  exhibit  predominantly  1  or  2  crossings. 


Table  3.  13.9  kcal/mol  Barrier  Crossings  for  Cr  +  CH3CI  in  Water 

Number  of  Crossings 

Initial  and  Final  States  of  Trajectories 

Reactant  •  Reactant 

Reactant  -  Product 

Product  -  Reactant 

Product  -  Product 

1 

0 

148 

0 

0 

2 

130 

0 

0 

110 

3 

0 

4 

6 

0 

4 

0 

0 

0 

2 

Of  the  63%  of  trajectories  which  cross  the  barrier  more  than  once,  most  recross  only  a  single  time. 
There  is  also  a  very  small  fraction  of  trajectories  which  have  3  or  4  crossings  of  the  energy  barrier. 
(Similar  recrossing  patterns  are  observed  for  35  trajectories  in  the  larger  263  water  system  in  which 
60%  of  the  trajectories  recrossed  the  barrier  once,  and  no  trajectories  crossed  more  than  twice.)  This 
observation  of  significant  recrossing  implies  a)  that  there  is  a  serious  deviation  from  simple  Transition 
State  Theory  (TST),5-6-29  which  assumes  all  activated  reactions  proceed  toward  products  with  no 
recrossings  of  die  barrier,46  and  b)  that  the  accepted  fundamental  reaction  "mechanism"  along  an  aver¬ 
age  free  energy  reaction  profile5-6  needs  to  be  reexamined. 

The  correction  to  the  TST  rate  constant  iP7  due  to  recrossings  of  the  transition  barrier  is  quanti¬ 
tatively  given  by  the  transmission  coefficient 


K  = 


(3.1) 


where  k  is  the  true  rate  constant.  We  will  compute  only  the  correction  to  i737  rather  than  the  rate  con¬ 
stant  itself,  as  the  latter  requires  a  knowledge  of  the  free  energy  of  activation  which  is  a  computation¬ 
ally  intensive  calculation  to  perform  using  molecular  dynamics.32-47-48  In  addition,  it  is  k  that  is  a 
direct  measure  of  the  dynamical  influence  of  the  solvent  in  affecting  the  rate  constant.  Under  the  con¬ 
ditions  in  which  trajectories  are  initialized  on  the  transition  barrier  surface  with  initially  forward 
momentum  in  the  asymmetric  stretch  reaction  coordinate,  then  run  forward  and  backward  in  time  from 
the  initial  conditions  to  determine  the  reaction  beginning  and  outcome,  the  transmission  coefficient  k  is 
given  by  the  Stable  States  Picture  of  reaction  dynamics15-49-50  as 

k  {  j*  «W01  )r~(U  9[x(-f)]  )j( 

Kmv * - nr. - •  a2) 


where  j+  represents  die  initially  positive  flux  across  the  transition  barrier  at  time  zero,  0[x(r)]  is  a  step 
function  for  forward  time  dynamics  which  equals  one  on  the  product  side  and  zero  on  the  reactant  side, 
x(t)  is  the  asymmetric  stretch  reaction  coordinate  defined  to  be  zero  at  the  transition  barrier,  8[x(-/)]  is 
a  step  function  for  reverse  time  dynamics  which  equals  one  on  the  product  side  and  zero  on  the  reactant 
side,  and  finally  (  )*  denotes  an  average  over  the  equilibrium  distribution,  normalized  by  the  reactant 
distribution.  (See  Ref.  (15)  for  a  more  detailed  discussion.)  Each  trajectory  i  in  the  equilibrium  distri¬ 
bution  has  a  particular  probability  w,  and  some  initial  velocity  v;  along  the  asymmetric  reaction  coordi¬ 
nate  and  is  integrated  out  to  0.5  ps  to  determine  the  reaction  fate.  Thus,  for  N  trajectories  in  the  ensem- 
;  ble,  k  translates  to,  and  is  computed  as15 

I  " 

Z  "i  |V,  |  Qi 

k  =  ~  ,  (3.3) 

i 
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where  +  indicates  all  trajectories  have  initially  positive  velocity  across  the  transition  state  surface  and 
the  value  of  Q,  depends  on  the  initial  and  final  states  of  the  trajectory  », 


+1 

if  Reactant  — ►  Product 

0 

-1 

if  Reactant  — *  Reactant  or  Product  — »  Product 

if  Product  ->  Reactant . 

(3.4) 

For  the  400  trajectories  in  the  ensemble  we  compute  k  to  be  0 .55.  This  represents  a  significant 
dynamical  correction  to  the  TST  rate  constant. 

2.  Recrossing  Pattern  Analysis 

Figure  5  qualitatively  illustrates  the  three  predominant  types  of  crossings  observed:  a)  direct  suc¬ 
cessful  reactant  -»  product  (RP)  transition  with  no  recrossing  (this  is  the  only  trajectory  type  in  the  TST 
description);  b)  a  single  recrossing  reactant  -*  product  -»  reactant  (RR)  after  the  transition  state  is 
crossed  in  the  forward  direction  and  c)  a  single  recrossing  product  -»  reactant  — >  product  (PP),  in 
which  the  recrossing  occurs  prior  to  an  initial  reactant  -»  product  crossing.  Why  are  there  recrossings, 
and  what  is  different  for  1  crossing  vs  2  crossing  cases?  We  find  that  these  features  can  be  correlated 
with  initial  solvent  configurations  at  the  transition  barrier  and  with  the  strength  of  the  solvent-solute 
interaction. 

A  first  indication  of  the  strength  of  the  reaction  system-solvent  coupling  is  revealed  by  examina¬ 
tion  of  the  time  dependent  friction,11'13’50  £(<),  which  the  solvent  exerts  on  the  reaction  coordinate  in 
the  neighborhood  of  the  transition  barrier.  This  is  computed  by  fixing  the  reaction  system  coordinates  in 
a  linear  configuration  at  the  gas-phase  transition  barrier  and  computing  the  time  correlation  of  the 
fluctuating  solvent  forces  resolved  onto  the  normal-mode  asymmetric-stretch  reaction  coordinate.15’51 
More  explicitly,  the  time  dependent  friction  is  defined  as 

CO)  =  (P/d)  <  FF(i)  )  ,  (3.5) 

in  which  F  is  the  force  acting  on  the  asymmetric-stretch  reaction  coordinate  with  effective  mass  p  and 
P"1  s  kBT.  Figure  6  shows  that  this  force-force  time  correlation  function  has  an  initially  fast  decay  out 
to  0.1  ps,  followed  by  a  short  plateau  region  of  duration  approximately  0.025  ps,  and  then  a  slow  decay 
toward  zero.  When  the  time  dependent  friction  is  computed  on  a  typical  bent  geometry  of  10°  at  the 
transition  barrier,  the  result  is  the  same  as  for  the  linear-geometry  friction.  We  will  return  to  this  time 
dependence,  as  associated  with  solvent  dynamics,  but  the  major  point  here  is  that  the  initial  magnitude 
of  the  friction 

C(r  =  0)  =  (p/p)  (  F2)  *  a>\  (3.6) 

indicates  a  very  strong  asymmetric  stretch-solvent  coupling,  u>;  =  890  cm-1,  to  be  compared  to  the  bar¬ 
rier  frequency  for  the  gas-phase  potential  tot  of  460  cm-1.  Thus,  it  is  not  surprising  that  this  strongly 
coupled  solvent  can  induce  recrossings  even  for  a  sharp  intrinsic  barrier.11'13’50 

We  now  investigate  the  source  of  this  strong  coupling.  A  comparison  calculation  of  c>;  for  an 
isolated  Cl"  ion  in  water  gives  0);  =  180  cm"1,  so  it  seems  reasonable  to  associate  the  strong  solvent 
coupling  noted  above  with  a  feature  of  the  reaction  system.  The  most  likely  association  is  with  the 
charge  migration  in  the  reaction  system  along  the  reaction  coordinate  (Fig.  2).  To  verify  this,  a  calcula¬ 
tion  of  the  time  dependent  friction  on  the  reaction  coordinate  is  made  in  which  the  charge  migration 
along  the  reaction  coordinate  is  made  to  pass  through  an  inflection  point  (cf  Fig.  7)  at  the  transition  bar¬ 
rier,  thus  giving  a  partial  charge  configuration.  Cl  *" — CHj6* — Cl  *"  which  is  the  same  as  before  but 
which  now  has  aero  first  derivative  with  respect  to  the  reaction  coordinate  at  the  transition  barrier.  This 
calculation  gives  0);  as  220  cm"1,  a  significantly  smaller  value  compared  to  the  890  cm"1  value  above, 
n  tnus  appears  that  the  rate  of  change  of  charge  along  the  reaction  coordinate  is  indeed  the  key  to  the 
strong  solvent  forces  experienced  by  the  reaction  system.  We  anticipated  this  result  in  Sec.  D,  where  it 
was  explained  that  an  extra  force  between  the  reaction  system  and  the  solvent,  in  addition  to  the  regular 
Coulombic  force,  arises  from  the  migration  of  charge  with  motion  along  the  reaction  coordinate. 
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The  critical  character  of  the  shifting  charge  distribution  is  further  revealed  by  correlating  the  reac¬ 
tion  fate  with  the  Coulombic  energy  of  interaction  of  the  two  chlorine  atoms  with  the  surrounding  sol¬ 
vent  at  the  transition  barrier.  For  the  400  trajectories  calculated,  we  find  that  the  reaction  outcome  is 
laigely  determined  by  which  of  the  two  chlorine  atoms  is  best  solvated  by  the  particular  configuration 
of  the  surrounding  solvent  molecules  at  the  transition  barrier.  For  example,  if  the  'left-hand”  side 
chlorine  has  a  more  favorable  interaction  energy  with  the  water  as  compared  to  the  'right-hand"  side 
chlorine,  then  the  reaction  tends  to  favor  forming  a  "left-hand”  side  chloride  ion,  and  this  is  true  regard¬ 
less  of  which  way  the  momentum  is  initially  headed  along  the  reaction  coordinate.  The  same  is  true  for 
the  opposite  situation,  in  which  the  "right-hand”  chlorine  atom  has  some  biased  attractive  potential 
energy  with  the  solvent  to  form  the  ion,  and  does.  In  these  two  cases  the  trajectory  recrosses  the  transi¬ 
tion  barrier  once  and  forms  the  ion  which  is  best  solvated  by  the  current  solvent  configuration.  If  pro¬ 
ducts  are  favored  then  products  are  formed,  and  the  product  side  is  where  the  reaction  begins  (cf  Fig 
Sc).  If  reactants  are  favored,  then  reactants  are  formed,  and  reactants  are  also  where  the  reaction  begins 
(cf  Fig.  Sb).  Finally,  we  find  that,  if  there  is  a  sufficiently  symmetric  solvent  potential  with  respect  to 
either  chlorine  atom  at  die  transition  barrier,  as  if  there  were  no  particular  preference  by  the  solvent  to 
have  either  the  "left"  or  "right-hand"  chloride  ion  formed  then  the  reaction  almost  always  proceeds  from 
reactants  to  formation  of  products  with  no  recrossing  of  the  transition  barrier  (cf  Fig.  Sa). 

Figure  8  illustrates  these  points  in  detail  by  showing,  for  ensembles  of  each  of  the  three  reaction 
outcomes  described  above  and  in  Fig.  S,  die  average  electric  potential  at  the  transition  barrier  at  the 
"left-hand"  chlorine  atom.  Cl,  compared  to  the  "right-hand"  chlorine  atom.  Cl',  as  a  function  of  the 
radius  of  neighbor  atoms  included  in  the  potential  calculation.  In  addition,  the  average  electric  poten¬ 
tial  is  computed  for  the  sum  of  all  three  ensembles.  This  electric  potential  V  is  computed  for  a  given 
radius  r  and  a  given  chlorine  atom  i  as 

V,<r)=  I  f  .  (3.7) 

i'a  s  ')  r<i 

where  is  the  distance  between  die  chlorine  atom  i  and  solvent  atom  j  and  is  the  charge  on  solvent 
atom  j.  Reactants  are  defined  as  Cl  +  CH3CT,  and  products  are  defined  as  CICHj  +  Cl'.  Figure  8 
shows  that  (a)  there  is  a  greater  electric  potential  at  "right-hand”  chlorine  atom.  O',  for  those  reactions 
which  cross  the  transition-state  barrier  twice  from  products  to  products  (PP),  and  (b)  a  larger  potential 
at  G  for  reactions  crossing  the  transition-state  barrier  twice  from  reactants  to  reactants  (RR),  but  (c)  no 
obvious  difference  in  electric  potential  for  die  reactions  which  cross  once  from  reactants  to  products 
(RP).  Note  that  the  average  for  all  three  reaction  cases  shows  no  average  bias  in  favor  of  one  of  the 
chlorine  atoms  at  die  transition  barrier,  as  is  appropriate  for  an  equilibrium  ensemble. 

Figure  8  also  shows  that  the  largest  electrical  potential  difference  for  the  PP  and  RR  reactions 
occurs  between  2.S  and  3.5  A  radially  outward  from  the  chlorine  atoms  and  subsides  to  a  nearly  con¬ 
stant  value  thereafter.  The  onset  at  small  neighbor  distance  of  differences  between  the  PP,  RP,  and  RR 
cases  shows  that  the  local  solvent  configuration  about  the  chloride  ions  largely  determines  the  reaction 
outcome.  A  further  indication  of  this  is  that  the  reaction  dynamics  are  essentially  the  same  for  the  small 
solvent  system  of  64  waters  and  the  large  263  water  system;  the  nearby  electrical  potentials  are  very 
similar  for  these  two  systems,  and  the  fairly  constant  potentials  at  longer  distance  from  the  transition 
barrier  apparently  have  little  effect  on  the  reaction  fate. 

All  this  suggests  that  the  transmission  coefficient  can  be  approximately  and  usefully  written  as  the 
contributions  arising  from  one  and  two  crossings, 

k  =  y-rr-  \(  u  0WO)  )Jr}  +  ( i*  eWOl  )JT")  -  < k  ®W-0]  >^>1 

.  ( j*  8[«(Q]  )jr° 

<>♦>* 

=  Prob(jym) ,  (3.8) 

where  (sym)  schematically  indicates  symmetric  solvent  configurations  which  always  give  rise  to  a 
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reactive  trajectory,  i*.  1  crossing,  while  (asym)  indicates  asymmetric  solvent  configurations  giving  rise 
to  2  crossings  and  which  never  result  in  an  overall  reaction  and  thus  yield  zero  net  contribution  to  k. 
This  reduces  k  to  the  first  term  which  we  have  approximately  expressed  as  Prob(jy/n),  the  relative  pro¬ 
bability  for  the  ’sufficiently”  symmetric  solvent  configurations.  Since  that  probability  is  less  than  unity, 
so  is  K. 

3.  Frozen  Solvent  Nonadiabatic  Solvation  Model 

Ultimately,  it  is  the  force  exerted  by  the  solvent  on  the  reaction  coordinate  plus  die  internal  force 
of  the  reaction  system  itself  which  drives  the  mat  Oon  toward  its  destiny.  We  have  found  above  that  the 
solvent  potential  energy  bias  between  the  two  chlorine  atoms  at  the  transition  barrier  is  a  strong  predic¬ 
tor  of  die  outcome  of  the  reaction.  This  strongly  suggests  that  there  is  a  force  F,  due  to  biased  solvent 
configurations,  which  may  shift  the  location  of  the  actual  transition  barrier  and  which  acts  along  the 
reaction  coordinate  to  drive  the  reaction  one  way  or  another.  Little  explicit  information  on  a  particular 
solvent  configuration  is  obtained  for  this  biased  force  from  the  calculation  of  the  solvent  friction  using 
Eq.  (3.5),  since  this  gives  an  average  square  force  for  an  ensemble  of  initial  conditions.  We  can,  how¬ 
ever,  examine  the  nature  of  the  ’true"  total  potential  which  defines  the  actual  transition  barrier  for  a 
particular  solvent  configuration  as  a  perturbation  from  the  gas-phase  transition  barrier.  While  a  more 
detailed  analysis  is  reserved  for  a  companion  paper,22  a  simple  analysis  suffices  here  to  expose  the  key 
features.  We  begin  by  assuming  a  nonadiabatic  frozen  solvent  model,  i.e.  that  the  solvent  does  not 
move  appreciably  on  the  time  scale  -0.02  ps  of  the  motion  which  decides  the  fate  of  the  reaction  pro¬ 
cess.  Next  the  "true"  potential  Vj(R)  as  a  function  ol  the  reaction  coordinate  R  which  contains  die 
potential  due  to  the  solvent  plus  that  of  the  gas-phase  reaction  system,  is  expanded  in  a  Taylor  series 
about  the  gas-phase  transition  barrier  location,  R  -  0, 

—  Rm  d"VT 

09> 

For  simplicity,  terms  with  order  greater  than  two  are  discarded  here.  The  zeroth  order  term  Vj<0)  is  die 
gas-phase  potential  energy  at  the  transition  barrier,  modified  by  the  interaction  with  the  equilibrated  sol¬ 
vent  which  lowers  that  energy  (Fig.  4).  Since  Vt<0)  is  already  incorporated  into  the  probability  that  we 
are  at  the  transition  barrier  and  that  the  molecular  dynamics  has  selected  the  particular  solvent 
configuration,  we  write  the  remaining  effective  potential  referenced  with  respect  to  it  as 

3Vr  i  ,  #VT 

<v«-*-ar  +  2  (3"» 

The  first  derivative  of  the  gas-phase  potential  at  R  -  0  is  zero  at  the  saddle  point.  Thus  only  the  frozen 
solvent  force  F,  contributes  to  the  first  derivative  term, 

dVT 

f'  =  ~!r  ’  (311) 

while  the  second  derivative  contains  contributions  both  from  the  gas-phase  potential  and  from  the 
frozen  solvent, 

1  - 1  , 

T1F=  2t“W  (312) 

in  which  is  a  nonadiabatic  frequency  describing  the  force  gradient  acting  on  the  reaction  coordi¬ 
nate  for  a  frozen  solvent  The  square  nonadiabatic  frequency  is  composed22  of  two  contributions 

“m.  =  w*  -  (Aw)2 .  (3.13) 

involving  the  bare  gas-phase  barrier  frequency  to*  and  a  solvent  bias  frequency  do.  The  latter  depends 
on  the  solvent  configuration  and  will  modify22  the  bare  barrier.  The  solvent  bias  frequency  reflects  a 
solvent  shift  in  the  barrier  curvature,  as  opposed  to  the  solvent  bias  force  F,  which  shifts  the  barrier 
position  (and,  as  we  will  see  presently,  its  energy).  As  shown  elsewhere,22  the  initial  value  of  the  time 
dependent  friction  is  related  to  the  average  square  solvent  bias  force 
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C(0)  =  (p/ji)  {  F* )  =  co£  ,  (3.14) 

while  (Ato)2  depends  on  instantaneous  solvent  configurations. 

The  effective  potential  Eq.  (3.10)  at  the  gas-phase  transition  barrier  can  now  be  rewritten  in  the 
suggestive  form 


V*>  =  --JM* a&P-W 


1  2 
-  -  2 


R  + 


F, 


wfij.  J 


F) 


(3.15) 


This  shows  tiiat  the  location  R  of  the  instantaneous  solvent  dependent  transition  barrier  is  shifted  away 
from  the  equilibrium  transition  barrier  (which  corresponds  in  location  to  that  of  the  gas-phase  transition 
barrier  if  ■  0)  by  an  amount 


Aif  =  - 


F, 


(3.16) 


while  the  shift  in  energy  from  the  equilibrium  transition  barrier  value  V^O)  is 
F] 


AV  = 


(3.17) 


These  features  are  illustrated  in  Fig.  9,  in  which  the  solvent  bias  potential  shifts  the  barrier  by  A/f  in 
position  along  the  reaction  coordinate  and,  very  importantly,  by  AV  in  energy  relative  to  the  equilibrium 
transition  barrier.  The  shift  in  energy  is  always  greater  than  or  equal  to  zero  with  respect  to  the  equili¬ 
brium  transition  barrier,  while  the  transition  barrier  will  shift  left  or  right  depending  on  the  nature  of  the 
particular  solvent  configuration.  We  will  presently  employ  this  simple  description  to  characterize  both 
the  crossing  patterns  and  k  in  detail. 

There  is  a  central  assumption  in  our  frozen  solvent  description  above  that  the  solvent 
configuration  does  not  undergo  any  appreciable  change  in  the  time  required  for  the  reaction  fate  to  be 
decided.  This  is  in  fact  reasonable,  considering  that  the  characteristic  time  scale50  on  a  barrier  of  fre¬ 
quency  (■)*  =  461  cm-1  is  2o)i'  »  0.02  ps.  (2 a>il  is  the  time  required  for  a  thermal  trajectory  to  decrease 
its  potential  energy  2 ktT  below  the  barrier  top.)  On  this  short  time  scale,  the  solvent  does  not  move  into 
a  significantly  different  configuration,  and  any  bias  (or  lack  of  bias)  in  potential  with  respect  to  the  two 
chlorine  atoms  does  not  break  down  appreciably  before  the  reaction  fate  is  decided.  This  is  also  in 
accord  with  the  minor  variation  of  £(r)  on  a  time  scale  of  approximately  1(T2  ps  (c/Fig.  6).  As  this  bias 
is  primarily  associated  with  the  charge-dipole  electric  potential  ( cf  Fig.  8),  the  only  thing  which  could 
break  the  bias  would  be  a  reorientation  or  a  significant  change  in  the  water  dipoles  within  the  time 
required  for  barrier  crossing.  The  only  solvent  motions  which  can  occur  on  this  brief  time  scale  are  the 
internal  vibrations  described  by  the  symmetric  and  asymmetric  stretches  and  bending  motions  of  water, 
with  normal  mode  frequencies  of  3832  cm-1,  3943  cm-1  and  1649  cm-1  respectively.52  However,  these 
vibrations  do  not  affect  the  orientations  of  the  water  dipoles  as  much  as  water  molecule  rotational 
motions  do,  and  the  gross  rotational  motions  are  slow  on  the  time  scale  for  barrier  crossing. 

Both  of  these  aspects  are  illustrated  in  the  time  dependent  friction  Fig.  6  and  in  Fig.  10,  which 
(hows  a  calculation  of  the  rotational  time  correlation  functions  for  pure  flexible  water  solvent  at  298  K. 
A  comparison  of  Fig.  10  with  the  time  dependent  friction  correlation  function  in  Fig.  6  shows  a  remark¬ 
able  similarity  between  the  two  characteristic  decay  times  and  profiles.  Both  correlation  functions  have 
an  initial  fast  decay,  followed  by  a  small  rise,  and  then  a  long,  slow  decay.  It  would  appear  likely  that 
the  major  time  development  of  the  time  dependent  friction  is  strongly  related  to  the  rotational  motions 
of  the  water  solvent,  and  this  seems  reasonable  in  light  of  the  observatK  ,i  that  the  strong  coupling 
between  the  solvent  and  reactions  system  is  largely  due  to  the  charge-dipole  forces,  which  change  most 
significantly  with  water  rotational  motions.  Thus,  the  rotational  motions  of  the  waters  do  not  effect  a 
major  change  in  the  solvent  configuration  on  the  time  scale  for  barrier  crossing;  rather,  the  solvent  is 
seen  to  move  on  a  much  longer  time  scale  after  the  decision  has  been  made  concerning  the  reaction 
outcome  by  the  instantaneous  configurations  of  the  water  molecules. 
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4.  Analysis  of  x 

The  features  of  the  solvent  bias  described  above  are  revealed  in  the  dynamics  because  those  tra¬ 
jectories  which  have  an  asymmetric  solvent  configuration  and  recross  the  gas-phase  transition  state  bar¬ 
rier  are  seeing  an  effective  potential  in  which  the  reaction  system  is  not  really  at  the  "true”  instantane¬ 
ous  solvent  configuration  dependent  transition  barrier  defined  by  the  effective  potential  [Eq.  (3.15)]: 
rather  they  are  shifted  from  that  location  by  an  amount  A#  in  position  [Eq.  (3.16)]  and  different  in 
energy  by  an  amount  AV  (Eq.  (3.17)].  Thus,  the  particular  trajectory  must  either  have  extra  kinetic 
energy  along  the  reaction  coordinate  equal  to  AV  in  order  to  react,  or  face  an  effective  barrier  which 
sends  it  back  («'.«.,  recrosses)  down  the  potential.  On  the  other  hand,  those  trajectories  which  have  a  sol¬ 
vent  configuration  which  produces  a  solvent  potential  which  does  not  lead  to  an  appreciable  AV  will 
usually  have  enough  kinetic  energy  to  react  without  recrossing.  This  suggests  the  following  analysis. 

Since  the  effective  potential  is  shifted  higher  in  energy  by  an  amount  AV  depending  on  solvent 
configuration,  a  trajectory  would  require  at  least  that  much  kinetic  energy  in  the  asymmetric  stretch  to 
overcome  the  barrier  and  make  a  successful  crossing.  This  suggests  that  a  plot  of  the  reaction  outcome 
for  a  given  kinetic  energy  in  the  reaction  coordinate  vs  AV  should  show  that  those  trajectories  which 
have  kinetic  energy  in  excess  of  AV  will  be  successful,  while  those  with  less  kinetic  energy  will  not 
This  idea  is  investigated  in  a  calculation  in  which  a  subset  of  100  initial  reaction  system  configurations 
are  randomly  selected  from  the  ensemble  of  400  initial  conditions  and,  for  each  initial  condition, 
random  positive  values  of  the  asymmetric  stretch  momentum  are  chosen  from  a  uniform  energy  distri¬ 
bution.  The  trajectories  are  then  run  long  enough  to  determine  whether  or  not  they  will  recross  the  bar¬ 
rier.  These  1000  trajectories  are  represented  as  points  on  a  "correlation”  plot  in  Fig.  11,  which  shows 
the  reaction  outcome  for  the  various  values  of  kinetic  energy  in  the  reaction  coordinate  vs  the  additional 
frozen  solvent  corrected  barrer  AV.  Here,  individual  trajectories  are  labeled  with  ’+’,  'O',  or  'A', 
depending  on  whetrherthey  cross  the  transition  barrier  once  from  reactants  to  products  (RP),  cross  twice 
from  reactants  to  reactants  (RR,  A R  >  0),  or  cross  twice  from  products  to  products  (PP,  A R  <  0),  respec¬ 
tively.  These  different  cases,  according  to  the  frozen  solvent  model,  should  be,  and  mostly  are  parti¬ 
tioned  into  regions  which  are  separated  by  a  line  of  slope  1,  dividing  the  RP  and  RR  trajectories,  and  a 
line  of  slope  -1,  dividing  the  RP  and  PP  trajectories.  The  lines  of  slope  1  and  -1  indicate  the  threshold 
kinetic  energy  in  the  asymmetric  stretch  required  by  a  trajectory  to  cross  the  additional  barrier  of  height 
AV.  This  figure  shows  that  the  outcome  of  the  reaction  is  indeed  largely  determined  by  the  conditions 
of  the  reaction  system  and  the  instantaneous  configuration  of  the  solvent  at  the  gas-phase  transition  bar¬ 
rier.  Thus,  we  only  need  a  knowledge  of  the  initial  kinetic  energy  in  the  asymmetric  stretch  reaction 
coordinate  and  the  shifted  barrier  height  AV  (determined  by  the  gas-phase  potential  added  to  the  frozen 
solvent  potential  due  to  the  solvent  configuration  at  the  gas-phase  transition  barrier)  to  predict  what  the 
likely  outcome  of  the  reaction  will  be  without  the  need  to  run  dynamics  away  from  the  gas-phase  tran¬ 
sition  barrier. 

The  transmission  coefficient  k  can  also  be  predicted  without  running  any  dynamics,  except  on  the 
transition  barrier,  from  a  knowledge  of  the  ensemble  of  solvent  configurations  of  the  system  with  the 
reacting  atoms  located  on  the  transition  barrier.  This  was  first  indicated  by  Eq.  (3.8),  in  which  k  is 
related  to  the  relative  probability  for  finding  a  symmetric  solvent  configuration  about  the  two  chlorine 
atoms  at  the  transition  barrier.  Now  we  include  the  effect  of  having  kinetic  energy  in  the  asymmetric 
stretch  reaction  coordinate,  which  may  overcome  the  additional  barrier  height  AV  and  result  in  a  suc¬ 
cessful  reaction.  The  relevant  probability  for  a  successful  reaction  given  AV  is22 

P(reaction  |  AV)  =  e~^v  ,  (3.18) 

and  tc,  the  probability  for  a  successful  reaction  in  the  ensemble  of  initial  conditions  on  top  of  the  equili¬ 
brium  barrier,  is  the  average 

K  *=  P(reactkm)  =  <  e P(AV) )  ,  (3.19) 

in  which  P(AV)  is  the  probability  for  having  a  particular  biased  solvent  configuration  which  gives  rise  to 
an  energy  shift  AV.  This  average  can  be  performed  using  the  MD  ensemble  of  initial  conditions,  and 
using  Eqs.  (3.11),  (3.12)  and  (3.17)  to  determine  AV.  The  result  is  k  -  0.47  ±0.05  which  is  in  satisfac¬ 
tory  agreement  with  the  value  k  -  0.55  ±0.05  determined  by  direct  MD  simulation. 
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In  fact,  an  analytic  theory  for  K  in  the  frozen,  non  adiabatic  solvent  limit  presented  in  the  compan¬ 
ion  paper22  based  on  the  general  van  der  Zwan-Hynes  nonequilibrium  solvation  rate  theory12’ 13  predicts 

that 


1 


2p(AV~’  (3‘20) 

where  the  average  barrier  shift  is  {  AV ).  The  average  |S  {  AV )  is  determined  from  MD  simulation  at 
the  transition  barrier  to  be  1.92,  yielding  tc  -  0.45,  again  in  good  agreement. 

In  summary,  the  frozen  solvent,  nonadiabatic  perspective  appears  to  give  both  a  qualitatively  and 
quantitatively  accurate  picture  of  the  model  S»2  reaction  mechanism  and  rate. 

Finally,  we  have  stressed  throughout  this  paper  the  critical  importance  of  the  short  time  scale 
dynamics  for  the  reaction;  Eq.  (3.20)  for  k  is  an  explicit  example  of  this,  since  it  involves  only  a  static 
average  for  a  frozen  solvent.  It  is  of  interest  to  ask  what  the  prediction  of  Kramers  theory53  for  k  would 
be.  In  this  theory,  the  fully  time  integrated  friction  of  Fig.  6  enters  as  the  friction  constant 


(3-21) 


C  =  |  dt  p/p  (  FF(t) )  =  <£>\x, 

in  the  Kramers  transmission  coefficient, 
k  =  Vl  +  (^2(ll4^)1-  (^2©^,) 

=  H  .  (3.22) 

for  high  friction.  Here  ©j,^  is  the  mean  barrier  frequency22  which  we  estimate  as  roughly  990  cm-1. 
(This  value  is  from  our  companion  paper.22  )  The  friction  constant  is  approximately  £  »  to\  t,  while  ©j 
is  the  square  root  of  the  initial  friction  value  of  890  cm*1  from  Eq.  (3.6),  and  t  is  the  lifetime  of  the 
fence  correlations,  estimated  from  Fig.  6  as  t  =  1  ps.  This  gives  the  Kramers  prediction  of  k  =  0.007, 
catastrophically  far  below  the  actual  K  «  0.SS.  The  long  time  scale  integrated  friction  in  the  Kramers 
picture  is  able  to  induce  significant  amounts  of  barrier  recrossing,  resulting  in  a  diffusion  controlled  bar¬ 
rier  passage  and  a  minuscule  K  value.  Clearly,  the  long  time  scale  friction  is  completely  irrelevant  on 
the  short  time  scale  during  which  the  reaction  fate  is  decided.11'13-51-54 


5.  Energy  relaxation 

We  now  return  to  the  question  of  why  there  are  only  1  or  2  crossings  of  the  transition  barrier  (i.e. 
virtually  no  observed  instances  of  more  than  2  sequential  barrier  crossings).  Here  we  examine  the  time 
scales  for  which  the  incipient  products  are  energetically  stabilized  in  the  solvent  Figure  12  shows  the 
internal  energy  of  the  reaction  system  as  a  function  of  time,  in  which  the  energy  is  partitioned  into 
[Cl — CHj — Cl]"  potential  energy  plus  CH3 — Cl  vibrational  kinetic  energy,  CH3 — Cl  rotational  kinetic 
energy,  and  the  translational  kinetic  energies  of  the  CH3 — Cl  molecule  and  of  the  CI~  ion. 

The  key  feature  here  is  that  in  approximately  0.025  ps  before  and  after  the  transition  barrier,  the 
total  energy  in  the  reaction  system  increases  by  approximately  one  keal/mol,  and  then  is  rapidly  dissi¬ 
pated  into  the  solvent  at  somewhat  longer  times.  This  initial  increase  in  energy  and  subsequent  loss  is 
primarily  associated  with  the  CH3 — G  vibrational  kinetic  energy  (c/  Fig.  12).  Note  that  although  the 
figure  shows  the  [Cl — CH3 — Cl]*  potential  energy  plus  CH3 — Cl  vibrational  kinetic  energy,  the 
increase  in  energy  shortly  before  and  after  the  transition  barrier  is  associated  with  vibrational  kinetic 
energy,  and  not  due  to  any  increase  in  the  total  potential  energy.  It  is  another  example  of  the  effects  of 
the  shifting  charge  on  the  dynamics  as  the  incipient  Cl'  ion  accelerates  off  the  transition  barrier  into  a 
more  favorable  ionic  interaction  with  the  surrounding  waters,  initially  causing  the  CH3 — Cl  to  "form" 
more  rapidly  compared  to  the  gas  phase  and  thus  increasing  its  vibrational  kinetic  energy.  The  subse¬ 
quent  loss  of  excess  vibrational  energy  to  the  solvent  at  longer  times  (2  0.1  ps)  occurs  very  rapidly  (We 
have  verified  that  this  is  largely  independent  of  the  presence  of  the  charge  switching.).  This  is  in 
marked  contrast  to  the  energy  decay  patterns  observed  in  the  gas  phase  in  Fig.  3,  and  is  also  quite 
different  from  that  found  for  the  A  +  BC  reaction  in  rare  gas  solvents.15  Thus,  the  picture  is  that  the 
initial  solvent  configuration  can  result  in  no,  or  a  very  small  number  of  recrossings,  and  once  the 
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incipient  species  favored  by  the  extant  conditions  is  formed,  very  rapid  stabilization  ensues. 

The  observations  above  raise  interesting  questions  for  future  study  about  the  vibrational  relaxation 
of  polar  molecules  in  polar  solvents.  But  first  and  foremost  they  demonstrate  that  the  reaction  dynamics 
are  quintessendally  a  question  of  both  'intrinsic”  reaction  system  and  solvent  dynamics.  The  stabiliza¬ 
tion  observed  intimately  involves  die  solvent,  yet  the  time  scale  is  so  short  that  dynamics,  and  not 
equilibrium  solvation,  is  critical.  Our  results  show  that  the  relaxation  is  not  occurring  on  the  monotoni- 
cally  decreasing  mean  potential  ( cf  Fig.  4).  Rather,  rapid  stabilization  occurs  in  the  complex  well  on  a 
subpicosecond  time  scale  due  to  the  strong  interaction  with  the  solvent  It  is  only  subsequent  to  this  that 
solvent  dipole  reorientation  on  a  picosecond  time  scale  (cf  Fig.  10)  will  occur  to  ’solvate"  the  complex 
and  establish  the  mean  potential  features  of  Fig.  4  in  the  neighborhood  of  the  products. 

IV.  Variation  of  Solvent,  Barrier  Height  and  Charge  Switching 

In  this  section  we  explore  the  effect  of  variation  of  solvent,  barrier  height  and  charge  switching  to 
examine  the  reaction  dynamics  under  different  conditions  and  to  determine  a  range  of  reaction  condi¬ 
tions  over  which  the  nonadiabatic  frozen  solvent  model  applies,  and  some  conditions  for  which  this 
model  breaks  down.  Table  4  presents  the  results  of  this  section  (see  below)  and  of  Sec.  Ill  above. 

A.  Variation  of  Mass  of  Solvent  Hydrogens 

The  picture  of  this  SN2  reaction  which  emerges  is  that  the  reaction  dynamics,  insofar  as  the  reac¬ 
tion  outcome  is  concerned,  are  largely  independent  of  the  solvent  dynamics  for  all  but  the  least  sharp  (5 
kcal/mol)  barrier  (see  below).  We  can  test  this  if  we  slow  down  the  solvent  vibrational  motions  and 
rotational  reorientation  time  by  replacing  the  mass  of  the  hydrogens  on  water  by  the  mass  of  oxygen, 
but  leaving  all  other  variables  in  the  system  unchanged.  In  this  way  we  change  the  short  time  aspects  of 
the  time  dependent  friction,  i.e.  slow  down  the  initial  decay  in  the  correlation  function  (cf  Fig.  13),  but 
the  initial  value  of  the  friction  is  unchanged,  since  that  is  a  static  quality  independent  of  the  solvent 
mass.  Thus,  we  make  the  waters  even  more  "frozen"  by  slowing  down  their  rotational,  vibrational,  and 
translational  motions.  We  find  that  the  reaction  dynamics  (e.g.,  recrossing  patterns),  and  the  energy 
decay  patterns  remain  essentially  the  same  as  with  the  H20  case,  and  K  (cf  Table  4)  matches  the  frozen 
solvent,  nonadiabatic  model  prediction.  The  fact  that  our  reaction  dynamics  in  the  HzO  solvent  system 
behaves  so  similarly  to  the  reaction  dynamics  in  an  even  more  frozen  solvent  system  further  buttresses 
the  frozen  solvent  picture  we  have  developed. 

B.  Variation  of  Barrier  Height 

The  frozen  nonadiabatic  solvent  picture  of  the  reaction  suggests  two  predictions  for  the  conse¬ 
quences  of  varying  the  central  Sn2  barrier  height. 

First,  for  higher  and  sharper  barriers,  since  the  bare  barrier  frequency  (i>k  is  increased,  the  reaction 
time  scale  -2io£'  is  diminished,  and  the  frozen  solvent  scenario  should  be  even  more  accurate.  Further, 
the  stronger  intrinsic  driving  force  for  the  reaction  should  more  easily  dominate  solvent  opposition  and 
k  should  rise.  This  is  apparent  from  Eqs.  (3.13),  (3.17)  and  (3.20).  Increasing  to*  increases  the  nonadia¬ 
batic  frequency,  to| ^  ,  which  in  turn  decreases  the  average  energy  shift  {  AV)  (cf  Fig.  9).  Thus,  k  will 
increase  according  to  Eq.  (3.20).  This  prediction  is  confirmed  in  Table  4  for  a  31.9  kcal/mol  central  bar¬ 
rier. 

Second,  if  the  barrier  height  is  lowered,  concomitantly  lowering  the  bare  barrier  frequency,  we 
expect  k  to  diminish  in  the  nonadiabatic  solvent  picture.  But  unfortunately  the  entire  basis  for  the 
frozen  solvent  description  is  mooted  by  the  longer  time  scale  for  the  reaction  and  the  associated  motion 
of  the  solvent  molecules.  We  thus  expect  a  serious  breakdown  in  the  frozen  solvent  description,  and 
Table  4  amply  demonstrates  this  point  for  a  S  kcal/mol  barrier.  For  this  case  we  find  that  (cf  Eq. 
3.13)  is  in  fact  negative  for  a  significant  fraction  of  solvent  configurations.  This  means  that  there  is  a 
"polarization  cage"  along  the  reaction  coordinate  as  described  by  van  der  Zwan  and  Hynes,11'13  and  a 
different  analysis  is  required. 

We  intend  to  return  to  the  low  barrier  case  in  a  subsequent  study  and  analyze  the  role  of  the  sol¬ 
vent  motions  in  detail  for  the  case  in  which  they  are  on  the  same  time  scale  as  the  effective  motion 


along  the  reaction  coordinate. 


C.  Variation  of  the  Rate  of  Charge  Switching 

As  seen  in  Sec.  m,  die  variation  of  charge  along  the  reaction  coordinate  has  a  dominant  effect  on 
the  solvent  friction  and  on  the  reaction  dynamics.  In  a  further  series  of  calculations,  we  have  deter¬ 
mined  that  slowing  this  charge  variation  down  reduces  the  number  of  trajectories  which  recross,  the 
amplitude  of  the  time  dependent  friction,  and  the  energy  pickup  in  vibration  shortly  after  the  transition 
state. 

Here  we  discuss  explicitly  only  the  transmission  coefficient  k  as  a  function  of  the  rate  of  charge 
variation.  In  the  nonadiabatic  solvent  perspective,  k  is  determined  from  Eq.  (3.20)  by  the  average  sol¬ 
vent  induced  barrier  shift  {  AV )  which  by  Eq.  (3.17)  is  quadratically  sensitive  to  the  solvent  bias  force 
F,  .  For  reduced  rates  of  charge  variation,  this  force  will  diminish.  Thus,  (  AV )  will  decrease  and  k 
will  rise.  Similarly,  for  an  increased  rate  of  charge  switching,  k  is  expected  to  decrease. 

These  predictions  are  verified  in  Table  4.  The  rate  of  charge  switching  can  be  varied  by  adjust¬ 
ing  the  parameter  K]  in  die  charge  switching  function  (cf  Appendix  A).  Changing  this  parameter  by  a 
scale  factor  x  changes  the  rate  of  charge  switching  by  the  same  factor.  When  the  rate  of  switching  of 
the  charge  along  the  reaction  coordinate  is  decreased  by  a  factor  of  two,  K  is  0.73.  This  is  slightly 
larger  than  the  result  K  -  0.55  for  the  full  switching.  Clearly,  the  effect  is  fairly  robust  If  we  decrease 
the  charge  switching  rate  by  a  factor  of  eight  then  k  rises  to  0.87;  significantly  less  recrossing  occurs. 
A  correlation  plot  is  presented  for  this  case  in  Fig.  14.  This  diagram  shows  that  the  range  of  AV  for  this 
reduced  rate  of  charge  switching  is  much  less  than  that  for  the  normal  rate  of  charge  switching  correla¬ 
tion  (cf  Fig.  11),  as  predicted  by  die  above  argument  The  reaction  is  approaching  the  "weak  solvation 
limit”  >2. 13  Id  which  the  solute-solvent  interaction  is  sufficiendy  weak  so  that  passage  over  a  sharp  bar¬ 
rier  is  negligibly  perturbed.  (Calculations  for  the  case  of  no  charge  switching  close  to  the  transition 
state  give  similar  results  as  for  the  1/8’th  charge  switching,  tc  -  0.90.  This  calculation  requires  a 
different  functional  form  to  describe  die  charge  switching  (cf  Fig.  7).)  Again  the  analytic  prediction 
gives  an  excellent  description,  documented  in  Table  4. 

D.  Collapsed  Charge  Model 

All  of  our  reaction  dynamics  have  been  computed  for  the  united  atom  definition  in  which  the 
charges  on  the  hydrogen  atoms  and  the  carbon  atom  of  the  methyl  group  are  reduced  to  a  single  posi¬ 
tive  charge  on  the  central  atom.  In  forming  the  united  atom,  the  effect  of  the  hydrogen  atoms  extend¬ 
ing  into  the  surrounding  solvent  with  their  associated  charges  is  lost,  and  this  may  result  in  the  hydro¬ 
gens  of  the  water  solvent  moving  in  unusually  close  to  the  chlorine  atoms  at  the  transition  state  in  a 
effort  to  maximize  the  attractive  coulombic  energy.  Up  to  now,  all  dynamics  have  been  computed  for  a 
charge  distribution  which  has  a  slightly  lower  positive  charge  on  the  central  united  atom  than  would  be 
given  by  the  united  atom  approximation  to  Jorgensen’s  potentials8'9  (cf  Appendix  A  and  Fig.  2).  As 
explained  in  Appendix  A,  we  feel  that  this  is  a  reasonable  distribution  to  mimic  the  reaction  system- 
solvent  interaction  when  the  CH3  group  in  the  S*2  system  is  collapsed  to  a  united  atom.  (This  view  is 
supported  by  the  agreement  for  die  model  noted  in  Sec.  D.3  with  the  radial  distribution  functions  com¬ 
puted  by  Jorgensen.8'9  )  Still,  the  dynamics  without  this  slight  positive  charge  diminution  are  of 
interest  We  have  determined  that  the  dynamics  for  the  united  atom  system  without  this  slight  charge 
redirection  (cf  Fig.  IS)  are  similar  to  those  with  it  but  k  is  only  0.30.  In  addition,  there  are  "polariza¬ 
tion  cage"  effects1113  as  seen  for  the  low  5  kcal/mol  barrier  dynamics  (see  above),  and  the  frozen  sol¬ 
vent  picture  no  longer  applies.  This  further  reduction  in  k  is  expected,  given  the  increased  chlorine- 
solvent  interaction  likely  in  this  calculation  (cf  Appendix  A). 

V.  Conclusions 

Molecular  dynamics  have  been  computed  for  a  model  S\2  reaction  in  solution,  and  the  fundamen¬ 
tal  reaction  picture  is  found  to  be  significantly  different  from  that  painted  by  a  standard  free  energy 
equilibrium  solvation  picture  of  solvent  effects.  The  free  energy  description  treats  the  solvent  as  a  sys¬ 
tem  which  is  always  equilibrated  to  the  reaction  system  at  all  points  along  the  reaction  coordinate,  and 
the  solvent  accordingly  affects  only  the  free  energy  profile  of  the  reaction  coordinate  through 


Table  4.  Transmission  Coefficients  from  Molecular  Dynamics  and  Frozen  Solvent 


Barrier 

Height 

(kcal/mol) 

Barrier 

frequency 

(cm-1) 

Relative  Rate 
of  Charge 
Switching 

Solvent 

type  #  of  mol. 

Molecular  Dynamics 

K tin  #  of  traj.  error 

Frozen  Solvent  Model 

Km**  #  of  traj.  error 

13.9 

4612 

1.0 

H20 

64 

0.55 

400 

±0.05 

0.45 

400 

±0.05 

13.9 

4612 

1.0 

h2o 

263 

0.64 

50 

±0.13 

0.51 

50 

±0.13 

13.9 

4612 

1.0 

C^O 

64 

0.44 

100 

±0.10 

0.46 

100 

±0.08 

31.9 

893.6 

1.0 

H20 

64 

0.74 

100 

±0.09 

0.74 

100 

±0.09 

mm 

231.7 

1.0 

h2o 

64 

0.40 

100 

±0.10 

0** 

100 

— 

13.9 

4612 

2.0 

h2o 

64 

0.38 

400 

±0.05 

0.26 

400 

±0.04 

13.9 

4612 

0.5 

h2o 

64 

0.73 

400 

±0.04 

0.67 

400 

±0.06 

13.9 

4612 

0.125 

h2o 

64 

0.90 

400 

±0.03 

0.87 

400 

±0.03 

13.9 

4612 

0.0 

h2o 

64 

0.91 

400 

±0.03 

0.91 

400 

±0.02 

13.9 

461.2 

1.0*** 

h2o 

64 

0.31 

100 

±0.10 

0** 

100 

— 

•  Kb*  =  (l+RAV))-*,  Eq.  (3.20). 

**  A  serious  breakdown  of  the  theory,  see  text 
***  Collapsed  charge  model,  see  text 

The  estimated  errors  for  the  k*d’s  are  calculated  using  a  multinomial  cell  probability  analysis  and  a 
95%  confidence  interval,  the  Km*’s,  from  a  propagation-of-enors  method.55 


equilibrium  solvation  and  stabilization  of  the  intermediates.  It  also  implies  that  simple  TST  is  a  good 
description  for  the  reaction  since  only  the  free  energy  profile  of  the  reaction  coordinate  is  changed  by 
the  solvent 

We  find  serious  deviations  from  this  picture  in  our  model  study.  Transition  State  Theory  breaks 
down  markedly,  as  recrossings  of  the  transition  state  surface  are  observed.  The  reaction  outcome  is 
highly  dependent  on  the  local  configuration  of  the  solvent  at  the  transition  state;  in  fact,  it  is  in 
significant  measure  determined  by  that  configuration.  The  reaction  in  the  transition  barrier  neighborhood 
is  characterized  by  nonadiabadc  solvation,11'13  i.e.  frozen  solvent  configurations,  rather  than  the  adia¬ 
batic  equilibrium  solvation  envisaged  in  simple  TST. 

This  central  feature  has  two  sources.  The  first  is  that  the  time  scale  during  which  the  reaction  fate 
is  determined  is  significantly  shorter  than  the  characteristic  time  scale  of  solvent  reorientation.  Thus, 
reactions  occur  in  a  "frozen  solvent"  at  the  transition  barrier.  The  second  is  that  there  is  a  rapid  varia¬ 
tion  of  charge  distribution  along  the  reaction  coordinate  close  to  the  transition  barrier,  and  this  leads  to 
pronounced  coupling  between  the  SN2  system  and  the  water  solvent. 

The  location  of  the  transition  barrier  for  a  particular  solvent  configuration  can  shift  markedly  from 
the  gas  phase  barrier  in  height  and  position  along  the  reaction  coordinate.  This  shift  is  not  accounted 
for  in  the  equilibrium  solvation  free  energy  picture  which  assumes  that  the  transition  state  has  the  same 
location  as  in  the  gas  phase  reaction,  but  is  shifted  uniformly  down  in  energy. 
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Instead  we  find  that  there  are  solvent  configuration  dependent  transition  barriers  shifted  up  in 
energy  from  the  equilibrium  value.  These  shifts  have  been  correlated  with  asymmetric  "solvation”  in 
which  one  of  the  chlorines  in  the  S»2  transition  state  is  differentially  stabilized  by  the  solvent,  and  the 
connection  to  the  observed  recrossing  patterns  has  been  made  in  considerable  detail 

A  simple  frozen  solvent,  non  adiabatic  solvation  theory  for  predicting  the  outcome  of,  and  the 
transmission  coefficient  for  a  solution  reaction  based  on  the  solvent  configuration  at  fire  equilibrium 
transition  state  is  given,  in  which  there  is  no  need  to  compute  any  dynamics  except  at  the  transition  bar¬ 
rier.  This  theory  is  tested  over  a  wide  range  of  reaction  system  variables  and  found  to  be  in  good 
agreement  with  our  MD  simulation  results.  In  contrast,  a  Kramers  description  fails  dramatically.  We 
find,  however,  drat  this  simple  theory  breaks  down  (as  expected)  in  certain  regimes,  for  example,  with 
low  barrier  heights  and  barrier  frequencies.  The  full  analysis  of  these  regimes  is  a  subject  for  future 
investigation. 

We  find  in  our  model  system  that  vibrational  energy  transfer  of  the  nascent  products  to  the  polar 
solvent  is  quite  rapid.  This  is  another  example  that  on  the  molecular  level  the  mechanism  of  the  reac¬ 
tion  process  is  not  determined  by  motion  on  the  mean  potential  free  energy  surface,  but  rather  requires 
an  investigation  of  the  dynamics  to  reveal  its  essential  features. 

Finally,  our  predictions  of  a  transmission  coefficient  correction  k  =  0.SS  for  the  13.9  kcal/mol 
model  Cl"  +  CH3C1  reaction  in  H20  could  fairly  be  interpreted  as  having  only  a  modest  quantitative 
effect  on  the  overall  reaction  rate  constant, 561 57  whose  solvent  influence  is  dominated  by  the  large  sol¬ 
vation  free  energies  of  the  reactants,  products  and  transition  state.1'10  This  perfectly  valid  perspective, 
however,  does  not  take  into  account  that,  according  to  our  study,  the  molecular-level  "mechanism''  of 
the  reaction,  and  the  solvent  participation  therein  differs  fundamentally  from  the  standard  conception. 
Further,  it  is  likely  that  significantly  smaller  k  values  -  and  hence  larger  numerical  deviations  from  the 
equilibrium  solvation  TST  rate  constant  -  will  occur  for  lower  barrier  SN2  reactions  (cf  Sec.  IV.B).  We 
hope  to  report  in  the  future  on  studies  of  this  regime,  as  well  as  on  related  matters  such  as  alternate 
descriptions5*’59  of  the  charge  variation  along  the  reaction  coordinate. 
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APPENDIX  A 

The  non-bonded  interactions  between  the  water  solvent  and  the  solute,  [Cl - CH3 - Cl]-,  are 

described  by  a  Lennard-Jones  Coulomb  (LJC)  potential  in  which  the  charge  on  the  solvent  atoms  is 
fixed  while  the  charge  on  the  solute  atoms  varies  as  a  function  of  reaction  progress  to  model  the  charge 
switching  process.  This  total  non-bonded  potential  is 

V*  =  £  Vj®  .  (Al) 

>■1 


^  £U  »JJ  ' 


in  which  ■  r,  -  ty  is  the  distance  between  solute  atom  1  and  solvent  atom  j,  n  is  the  number  of  sol¬ 
vent  atoms,  the  three  solute  atoms  are  denoted  by  Cl  =  A,  CH3  =  B,  CT  =  C,  Qv  is  the  charge,  Ay  and 
£v  are  the  Lennard-Jones  constants,  and  VWau'  is  the  water  solvent  intermolecular  potential.21 

The  charge  and  Lennard-Jones  parameters,  Qr  are  composed  of  multiplicative  factors  of 

the  specific  solvent  and  solute  parameters  as  follows.  Let 

G?;  =  Mj .  (A3) 


Bq  =  bibj ,  (A5) 

when  q„  a„  b,  an  solute  parameters  and  qr  ar  bi  an  solvent  parameters,  then  Eq.  (A2)  becomes 

--^4+ . . . •r>>-,,r»+i •  •  ^ •  <A6) 

*«A  l  rV  *W  'I  J 

The  charge  switching  is  modeled  by  allowing  the  charge,  qb  on  the  solute  atoms  to  vary  as  a 
function  of  reaction  progress  in  a  fashion  similar  to  the  one-dimensional  charge  switching  prescription 
of  Jorgensen,*  although  our  reaction  progress  is  conveniently  defined  to  be  ('ab-'bc)  k>  allow  for  a  full 
three-dimensional  description  of  die  charge  switching.  In  addition  the  Lennard-Jones  parameters  a,  and 
can  also  be  made  variable  with  reaction  progress  such  that 

*  f.f'AB-'Bc)  •  (A7) 

■  a^fAB-ric) .  (A8) 


Bi  *  *,<'ab-,bc)  • 

The  force  on  the  solute  atom  due  to  solvent  interaction  is  given  by 
pNB  _  -dVNB 

*  ar* 
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_  V  +  b>  dbi } 

ar*  r?  dr*  "  ^  3r*  J 

i->  l  r*  rV  't  J  4 

>il  «»aL  rV  rv  M/  J 


Notice  that  the  solute  charge  and  Lennard-Jones  parameters  are  differentiated  with  respect  to  r*  since 
these  are  now  variable  with  the  reaction  progress.  These  terms  are  expanded  as  follows. 

For  g  =  qA,b\  ijt  =  AJB.C 

^r|c)  ^rAB_rBc)  =  8i’  ^t(rAB~rBc)  .  (All) 

in  which  the  explicit  gradients  of  (r^B-r|c)  are 

VA(dB-4c)  =  ^  ,  (A12) 

VB(dB-4c)  =  ^  .  (A13) 


^c(rAB~ric)  =  2tbc  • 
Definition  of  the  new  variables 
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(A16) 


allows  us  to  write  the  non-bond  forces  on  the  solute  atoms  due  to  the  yth  solvent  atom  in  terms  of  A j 
and  F„  as 

F5Jf  =  r^-  A/ab.  (A17) 

-  A/ca  »  (A18) 

and 


Fq  =  rqrQ  -  A/ac  • 


(A19) 


Hem,  Ay  is  of  the  same  fonn  as  the  first  term  of  except  that  qit  a„  and  £>,  are  replaced  by  2q',  2a', 
and  2b',  respectively.  Also  note  that  I"jjr#  is  just  the  normal  LJC  force.  Summing  over  all  the  solvent 
atom  interactions,  the  total  force  on  each  solute  atom  is  then 


-A/ab|. 

Fb®  =  £  Fg?  =  £  l^B/By  -  AyTcxj . 


(A20) 

(A21) 


and 


► 

f^  =  ZFq  =trc/Q 
>■> 


(A22) 


Our  adopted  functional  form  for  die  solute  LJC  parameters,  (<?„  a„  b,),  is,  for  g  ■  q#Jb  and  i 
gi  *  g.<ri b~4c) 

=  aif\(rAB~rBc)  +  $ifl(.rKB-rBc)  +  if 

=  aifi  +  Wi  +  i! . 
in  which  /,  is  an  odd  function 

Zi('AB-'ic)  =  tan'1K,(riB-»ic) . 
and/2  is  an  even  function 

,  ,  ■ 


A.B.C, 


(A23) 

(A24) 
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Similarly,  the  derivatives  far  g  =  i  =  A,B,C  are 


oJ/i'+Pft'. 


(A27) 


(A28) 


and 


(K2(rAB-ric)J  +  l)2 


(A29) 


We  choose  these  functions  rather  than  Jorgensen’s  functions8-9  because  of  the  non-differentiability  of 
die  latter  at  the  transition  state  as  well  as  the  requirement  for  defining  these  functions  for  the  three- 
dimensional  space  in  which  we  calculate  dynamics.  We  have  found  it  feasible  to  use  a  piecewise  fifth 
degree  polynomial  fit  to  these  functions  and  their  derivatives  when  we  actually  perform  the  dynamical 
calculations. 


V  -A-.- 


\  • 


The  central  atom  is  modeled  as  a  united  CH}  group,  which  is  given  the  diameter  of  a  methyl 
group  (4.25  A)  and  contains  no  explicit  hydrogens.  There  is  some  question  of  how  the  charges  on  the 
individual  hydrogen  and  carbon  atoms  should  be  summed  and  collapsed  down  to  form  a  single  partial 
positive  charge  on  the  central  united  atom.  In  forming  the  united  atom,  the  effect  of  having  the  hydro¬ 
gen  atoms  sticking  out  in  the  solvent  with  their  partial  positive  charges  is  lost,  and  this  may  allow  the 
hydrogen  atoms  of  the  water  molecules  to  move  unusually  close  to  the  two  chlorine  atoms  at  the  transi¬ 
tion  state  to  maximize  the  attractive  potential  energy.  To  compensate  for  the  loss  of  the  repulsion  force 
between  hydrogens  of  the  methyl  group  which  would  keep  the  hydrogens  of  the  waters  from  coming 
too  close  to  the  chlorine  atoms,  an  effective  reduction  in  the  negative  charge  of  the  chlorine  atoms  at 
the  transition  state  is  required  to  lower  the  attractive  potential.  Reducing  the  amount  of  negative  charge 
on  the  chlorine  atoms  at  the  transition  state  means  that  the  central  atom  must  carry  more  of  the  partial 
negative  charge  so  that  total  charge  conservation  is  maintained.  This  results  in  the  central  atom  having 
a  slightly  less  positive  charge  than  is  predicted  by  the  united  atom  definition  of  summing  the  charges  of 
the  methyl  and  collapsing  them  down  to  a  single  charge.  The  charge  switching  parameters  are  deter¬ 
mined  first  for  the  united  atom  approximation,  which  gives  an  overall  reasonable  fit  to  Jorgensen’s  one- 
dimensional  functions,*  with  K,  =  0.75  and  tc2  =  0.25.  Then  the  parameters  are  adjusted  so  that  the 
charge  on  the  central  group  decreases  (becomes  less  positive)  at  the  transition  state,  and  the  a?  parame¬ 
ters  are  chosen  to  match  the  dipole  moment  of  the  separated  CH3C1  molecule  as  described  by  Jorgen¬ 
sen.8  Figure  2  illustrates  the  variation  of  charge  on  the  three  atoms  of  the  reaction  system  as  a  function 
of  the  reaction  progress,  (rAB-rBc)*  Table  A1  gives  the  parameters.  The  charge  variation  for  the 
Jogensen  methyl  charges  collapsed  down  onto  a  single,  united  atom  are  shown  in  Fig.  15. 

The  nonvariable  charges  on  the  water  atoms  are  chosen  from  Watts,21  while  the  Lennard-Jones 
parameters,  A „  and  6U,  for  pure  solute-solute  and  pure  solvent-solvent  interactions,  are  taken  from 
Berendsen.60  The  combining  rules 

e„  =  (e*  Ci/)'4  *  (A30) 

and 

%  «  ~  (<*„  +  Of) ,  (A31) 

in  which  Au  *  4e„<r’J  and  Bu  =  are  used  to  compute  the  Lennard-Jones  parameters  for  solute- 

solvent  interactions.  For  convenience  in  our  calculations,  we  define  Ohh  for  hydrogen  atom  interactions 
to  be  equal  to  Oqo  for  oxygen  atom  interactions  (There  is  no  well  defined  parameter  for  H-H  and  this 
definition  is  used  for  consistency  in  the  combining  rules.).  Using  the  above  combining  rules,  there  is  a 
relationship  between  a„  b,  and  the  Lennard-Jones  parameters  for  pure  atom-atom  interactions  as  fol¬ 
lows: 

If  I  -  O.H  then 

a,  =  2£?o&  , 

b,  =  2eJ?Ooo  .  (A32) 

If  i  »  G,  CHj  then 

12 


(A33) 

The  parameters,  yf  and  yf  are  given  in  Table  Al.  To  focus  on  the  influence  of  the  switching  charges  on 
wic  dynamics,  die  constants  of,  JS“  a?  and  Jlf,  given  in  Table  Al,  are  set  equal  to  zero  to  turn  off  the 
variation  of  the  Lennard-Jones  parameters  with  reaction  progress. 


22  - 


Finally,  far  the  solvent  atoms,  the  forces  due  to  the  presence  of  the  other  solvent  molecules  and 
the  solute  atoms  are 


•T  *  £  r,r,  -  3V"‘ 


dry 


=  ryAr,A  +  ffiTjb  +  rycfyc  ~ 


dry 


(A34) 


Table  A1 


< 

;  /(kcal  A/ mol) 74 

a,  /(kcal  A!2/mo!)* 

b,  /(kcal  A^/moI)*4 

1 

a? 

P7 

a 

a?  p;  a 

P?  )? 

A 

-4.3688 

1.1255 

-11.3645 

0.0 

0.0 

6008.2526 

0.0 

0.0 

62.6221 

B 

0.0000 

-2.2510 

4.5020 

0.0 

0.0 

2623.0416 

0.0 

0.0 

47.1798 

C 

4.3688 

1.1255 

-11.3645 

0.0 

0.0 

6008.2526 

0.0 

0.0 

62.6221 

O 

0.0000 

0.0000 

-12.00666 

0.0 

0.0 

793.322 

0.0 

0.0 

25.0128 

H 

0.0000 

0.0000 

6.00333 

0.0 

0.0 

1245 

0.0 

0.0 

0.03935 

APPENDIX  B 

The  switching  function,  0(2),  satisfies  the  following  criteria: 
for  2^0  0(2)  =  1  o'(i)  =  0 


for  0  <  2  <  2q 
for  2  2  20 


0  <  0(2)  <  1  <j'(2)  <  0  . 
0(2)  =  0  o'  (2)  =  0 


(Bl) 


The  functional  form  of  0(2)  in  our  calculations  is 


0(2)  = 


1 


1  -  (2/20)3(10  -  (2/2o)(15  -  6(2/2o))) 
0 


tSO 

0  <  2  <  20  . 
a  >  *o 


(B2) 


Denoting  Q  by  A,  CHj  by  B,  and  Cl'  by  C,  and  letting  p  =  (r^a-rlc),  and  k  be  an  adjustable  positive 
constant,  the  image  points  of  the  complex  are 


r,  =  [l-o(Kp)]rA  +  o(Kp)rB 


(B3) 


and 


r2  =  a(-Kp)r„  +  [l-a(-Kp)]rc  .  (B4) 

For  p  »  0  we  have  A“  4  BC,  for  p  =  0  we  have  {A - B - C]‘,  and  for  p  «  0  we  have  AB  4  CT. 

Thus  we  have  a  smooth  transition  of  grouping  along  the  reaction  coordinate,  minimizing  the  effects  of 
our  finite  system  size. 

For  the  case  where  p  >  0,  the  modification  to  the  solute-solvent  potential  is  as  follows: 

Define 


image  point  of  group 


* 

a  = ' 

> 


ri 

f2 


a  =  A 
a  =  BC 


(B5) 


rp  ■  image  point  of  solvent  molecule  3  =  position  of  oxygen  atom  , 

tgp  =  ra  -  Tp  , 

12  ■  square  of  switching  region  minimum  value  , 


V„  *  Lennard-Jones  Coulomb  potential  between  solute  atom  i  and  solvent  atom  j 


(B6) 

(B7) 

(B8) 

(B9) 


The  potential  energy  of  interaction  between  solute  atom  i  and  solvent  atom  j,  where  atom  j  is  a  member 


> 


>1 

.,%1 


-ri 


If 

>  -1 


m 

.*> 


m 


'  *  *  •  '  *  j 


of  solvent  molecule  p,  is 

Va  =  CTlr^-j2)^  .  (BIO) 

The  seme  type  of  analysis  applies  when  p  <  0.  For  the  solvent-solvent  interactions,  the  interaction 

potential  between  two  atoms,  i  and  j,  from  different  solvent  molecules,  3)  and  3*.  is 

Vj,.  (Bll) 

For  our  system  of  64  waters,  i0  =  11.0  A2,  s  =  5.94  A,  and  k  =  0.9. 
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Figure  Captions 


Figure  1.  Gas  phase  potential  eneigy  LEPS  surface  contour  plot  as  a  function  of  the  AB  and  BC  bond 
lengths,  rt  and  /■*  for  a  linear  arrangement  of  [Cl — CH3 — Cl]-,  a  barrier  height  of  13.9  kcal/mol  with 
respect  to  the  bottom  of  die  wells,  and  a  barrier  frequency  of  a>b  -  461  cm-1.  The  contour  values  are  in 
kcal/mol  units,  and  the  zero  in  potential  is  far  separated  species  Cl-  +  CH3C1. 


Figure  2.  Variation  of  charge  on  the  three  atoms  of  the  reaction  system  as  a  function  of  -  r\c,  in 
which  the  charge  is  given  in  units  of  the  charge  of  an  electron.  Note  that  the  charge  variation  on  the 
central  united  atom,  CH3,  shows  a  slight  dip  (i.e.  becomes  less  positive)  in  charge  near  the  transition 
barrier.  This  is  different  than  the  slight  increase  in  positive  charge  which  would  be  obtained  if  one  sim¬ 
ply  summed  the  charges  on  the  atoms  of  CH3  to  a  single  point  charge,  as  is  discussed  in  Appendix  A 
and  shown  in  Fig.  IS. 


Figure  3.  Plot  of  the  average  internal  energy  redistribution  of  [Cl — CHj — Cl]",  versus  time  for  the 
reaction  Cl-  +  CH3C1  — >  CICHj  +  Cl-  on  the  13.9  kcal/mol  energy  barrier  surface  for  an  ensemble  of 
128  trajectories  in  the  gas  phase.  Time  zero  is  when  the  reaction  system  is  first  released  at  the  saddle 
point.  The  barrier  height  for  the  linear  configuration,  relative  to  separated  species,  is  approximately  3.6 
kcal/mol.  The  total  energy  is  partitioned  into  [Cl — CHj — CI]~  potential  energy,  [Cl — CH3— Cl]- 
potential  energy  plus  CH3— Cl  vibrational  kinetic  energy,  CH3 — Q  rotational  and  translational  kinetic 
energy,  and  the  translational  kinetic  energy  of  the  CC  ion.  The  vibrational  kinetic  energy  plus  total 
[Cl — CH3 — Cl]-  potential  energy  increases  to  a  maximum  in  about  0.25  ps  before  and  after  the  transi¬ 
tion  barrier  but  decreases  at  longer  times.  This  build  up  in  energy  and  subsequent  decrease  is  observed 
on  a  shorter  time  scale  for  the  translational  and  rotational  degrees  of  freedom.  The  kinetic  energy  at 
longer  times  is  seen  to  transform  into  potential  energy,  and  a  metastable  complex  is  formed.  In  the  gas 
phase,  in  contrast  to  the  solution  trajectories  described  later,  the  total  [Cl — CH3 — Cl]-  potential  plus 
kinetic  energy  is  conserved. 


Figure  4.  Schematic  diagram  showing  a  typical  thermodynamic  free  energy  profile  along  the  reaction 
coordinate  for  a  Sn2  reaction,  X-  +  RY  -»  XR  +  Y-,  in  the  gas  phase  (solid  line)  and  in  a  polar  solvent 
(dashed  line).  The  standard  free  energy  picture  is  that  the  polar  solvent  will  stabilize  the  higher  charge 
density  reactants  and  products  more  than  the  lower  charge  density  transition  state  complex.  The  energy 
wells  in  the  gas  phase  are  characteristic  of  a  metastable  chaige  dipole  complex  which  can  form  close  to 
the  transition  barrier,  and  will  mostly  disappear  in  solution  (in  thermodynamic  free  energy  terms). 


Figure  5.  Schematic  illustration  of  some  predominant  recrossing  patterns  observed  in  the  molecular 
dynamics  simulation  of  A  +  BC  reaction  dynamics  on  low  barriers  in  solution.  The  dividing  line 
represents  the  transition  barrier  dividing  surface  between  reactants  on  the  left  and  products  on  the  right, 
while  the  arrows  indicate  die  direction  of  the  trajectories:  (a)  is  a  direct,  successful  reactant  -»  product 
(RP)  transition  with  no  recrossing  (Only  this  type  of  trajectory  occurs  in  the  TST  picture.);  (b)  is  a  sin¬ 
gle  recrossing,  reactant  — »  product  — >  reactant  (RR),  after  the  transition  barrier  is  crossed  in  die  forward 
direction,  and  (c)  is  a  single  recrossing,  product  -»  reactant  -»  product  (PP),  in  which  the  recrossing 
occurs  prior  to  an  initial  reactant  -»  product  crossing. 


Figure  6.  Solvent  time  dependent  friction  on  the  reaction  coordinate  for  a  linear  geometry  of 
[Cl — CH3 — Gr  «  the  transition  barrier  as  function  of  time,  normalized  by  the  value  of  the  friction  at 


tune  zero.  The  magnitude  of  the  friction  at  time  zero  is  <o£,  with  o>;  =  890  cm-1.  When  the  time  depen¬ 
dent  friction  is  computed  for  a  typical  bent  geometry  of  [Cl — CH3 — Cl]-  at  the  transition  barrier,  die 
result  is  the  same  as  for  the  linear  geometry. 

Figure  7.  Variation  of  charge  on  the  three  atoms  of  the  reaction  system  as  a  function  of  rj*  -  in 
which  the  charge  is  given  in  units  of  the  charge  of  an  electron.  Here,  the  charge  migration  along  the 
reaction  coordinate  is  made  to  pass  through  an  inflection  point  at  the  transition  barrier,  thus  giving  a 
partial  charge  configuration  which  is  the  same  at  the  transition  barrier  as  in  Fig.  2  but  which  has  a  zero 
first  derivative  with  respect  to  the  reaction  coordinate  at  the  transition  barrier. 


Figure  8.  Illustration  of  the  coulombic  bias  at  the  transition  barrier  for  ensembles  of  each  of  the  three 
reaction  outcomes  described  in  Fig.  S.  The  average  electric  potential  (cf  Eq  3.7)  at  the  "left-hand" 
chloride  atom.  Cl,  is  compared  to  that  for  the  "right-hand"  chloride  atom,  CT,  as  a  function  of  the 
radius  of  a  sphere  of  neighboring  atoms  included  in  the  calculation  of  the  potential.  In  addition,  the 
average  electric  potential  is  computed  for  the  sum  of  all  three  ensembles,  i.e.  for  an  ensemble  of  essen¬ 
tially  all  trajectories.  Reactants  are  defined  as  Cl  +  CHjCT  and  products  are  defined  as  C1CH3  +  CT. 
The  figure  shows  that  there  is  a  greater  electric  potential  at  the  "left-hand”  chloride  atom.  Cl,  for  reac¬ 
tions  crossing  the  transition  barrier  twice  from  reactants  to  reactants  (RR),  and  a  larger  potential  at  die 
"right-hand"  chloride  atom.  Cl',  for  those  reactions  which  cross  the  transition  barrier  twice  from  pro¬ 
ducts  to  products  (PP),  but  no  obvious  difference  in  electric  potential  for  the  reactions  which  cross  once 
from  reactants  to  products  (RP),  nor  for  the  sum  of  all  trajectories. 


Figure  9.  Schematic  representation  of  the  effect  of  the  frozen  solvent  linear  and  quadratic  potentials 
when  added  to  the  gas  phase  potential  energy  profile  along  the  reaction  coordinate.  The  frozen  solvent 
potentials  shift  the  location  of  the  equilibrium  transition  barrier  forward  or  backward  along  the  reaction 
coordinate  and  to  higher  energy.  Note  that  the  peak  of  the  gas  phase  potential  is  lowered  to  V^O), 
which  includes  equilibrium  solvation. 


Figure  10.  Rotational  time  correlation  functions,  P,  and  Pj,  for  the  H-O-H  angle  bisector  for  87 
molecules  of  liquid  water.  These  correlation  functions  are  defined  as  (  Pm  [d-d(r)] ),  where  P„  is  the  n’th 
Legendre  polynomial,  and  d(<)  is  the  unit  vector  along  the  angle  bisector  (approximately  along  the 
dipole  moment  vector)  at  time  r.  Note  the  similarity  between  these  correlation  functions  and  the  time 
dependent  friction  for  water  acting  on  the  reaction  coordinate  as  shown  in  Fig.  6. 


Figure  11.  Outcomes  of  1000  trajectories  on  the  13.9  kcal/mol  barrier  and  normal  charge  switching  (c/ 
Fig.  2),  plotting  initial  kinetic  energy  along  the  asymmetric  stretch  reaction  coordinate  versus  the  rela¬ 
tive  shift  in  energy,  AV,  from  the  gas  phase  transition  barrier.  The  symbol  V  represents  trajectories 
which  cross  the  transition  barrier  once  from  reactants  to  products  (RP),  circles,  'O’,  represent  trajec¬ 
tories  which  begin  as  reactants  and  cross  the  transition  barrier  twice  to  end  up  as  reactants  (RR),  while 
triangles  ,'A',  are  trajectories  which  begin  as  products  and  cross  the  transition  barrier  twice  to  reform 
products  (PP).  The  lines  of  slope  1  and  -1  represent  the  threshold  of  kinetic  energy  required  to  over¬ 
come  the  additional  potential  energy,  AV,  to  give  a  successful  reaction. 


Figure  12.  Plot  of  die  average  energy  partitioning  of  [Cl — CH3 — Cl]-  versus  time  for  the  reaction 
Cl-  +  CHjCi  -»  ClCHj  +  CT  on  the  13.9  kcal/mol  energy  barrier  for  an  msemble  of  100  trajectories  in 
water  solvent.  Time  zero  is  when  the  reaction  system  is  first  released  he  saddle  point  The  total 
reaction  system  energy  is  partitioned  into  (Cl — CH3 — Cl]-  potential  energy  plus  CH3 — Cl  vibrational 
kinetic  energy,  CH3 — G  rotational  kinetic  energy,  and  the  translational  kinetic  energies  of  the  CH} — Cl 


and  of  the  Cl"  ion.  Note  that  the  increase  in  total  reaction  system  energy  shortly  before  and  after  the 
transition  barrier  is  largely  due  to  an  increase  in  the  kinetic  energy  of  vibration.  Note  that  the  vibra¬ 
tional  energy  of  CH3 — Cl  decays  rapidly  to  solvent  on  a  subpicosecond  time  scale. 


Figure  13.  Solvent  time  dependent  friction  on  the  reaction  coordinate  as  a  function  of  time,  normalized 
by  the  value  of  the  friction  at  time  zero,  for  a  heavy  mass  water  system  in  which  the  mass  of  hydrogens 
is  replaced  by  the  mass  of  oxygen,  but  all  other  variables  remain  unchanged.  The  magnitude  of  the  fric¬ 
tion  at  time  zero  is  the  same  as  for  the  HzO  solvent  system  (tty  =  890  cm'1),  but  the  initial  short  time 
decay  in  the  correlation  is  slower  than  that  for  pure  water  which  is  shown  in  Fig.  6.  It  is  interesting  to 
note  that  the  longer  time  behavior  is  nearly  the  same  as  in  Fig.  6. 


Figure  14.  Outcome  of  300  trajectories  on  the  13.9  kcal/mol  barrier  with  a  1/8’th  charge  switching  rate. 
Initial  kinetic  energy  along  the  asymmetric  stretch  reaction  coordinate  is  plotted  versus  the  relative  shift 
in  energy,  AV,  from  the  gas  phase  transition  barrier.  For  details,  see  Fig.  11. 


Figure  15.  Collapsed  charge  version  of  variation  of  charge  on  the  three  atoms  of  the  reaction  system  as 
a  function  of  r%g  -  r^,  in  which  the  charge  is  given  in  units  of  the  charge  of  an  electron.  The  slight 
increase  in  positive  charge  for  the  CH3  group  is  obtained  by  summing  the  charges  on  the  atoms  of  CH3 
to  a  single  point  charge,  as  is  discussed  in  Sec.  IV  and  in  Appendix  A. 
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